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1  Foreword 


Image  superresolution  refers  to  methods  that  increase  spatial  resolution  by  fusing  information 
from  a  sequence  of  images,  acquired  in  one  or  more  of  several  possible  ways.  The  high 
resolution  filtered  image  is  constructed  from  the  aliased  (undersampled)  noisy  and  blurred 
frames  with  either  subpixcl  shifts  or  the  use  of  intentional  blurring  by  designing  lenses  with 
different  point  spread  functions. 

A  demand  for  higher  resolution  is  seen  in  many  holds  including  bio-medical  imaging 
(for  purposes  like  image-guided  surgery  and  image-assisted  medical  diagnosis),  entertain¬ 
ment  (high  definition  television  or  HDTV),  satellite  and  astronomical  imaging,  chemical  and 
biological  research  (high  resolution  electron  microscopy),  military  surveillance  and  remote 
sensing.  Indeed,  this  demand,  in  many  cases,  exceeds  the  maximum  resolution  capability  of 
current  acquisition  systems.  In  other  words,  the  current  state  of  image  sensor  technology 
acts  as  a  limiting  factor  in  acquisition.  In  various  other  cases,  factors  such  as  cost,  physical 
attributes  like  size  and  weight,  and  quality,  instead  of  technology,  constrain  the  maximum 
resolution  obtainable  from  the  acquisition  device.  Yet  another  set  of  cases  exist  in  which 
resolution  is  compromised  due  to  the  need  for  flexibility  and  robustness  with  respect  to  var¬ 
ious  environmental  conditions,  among  others.  In  all  such  cases,  the  solution  to  the  need 
for  higher  resolution  necessitates  the  design  of  technological  devices  in  the  form  of  digital 
image  processing  algorithms  to  satisfy  the  demand  for  high  quality  and  high  resolution  (HR) 
images  and  video. 

2  Statement  of  the  Problem  Studied 

Wavelet  superresolution  is  a  topic  that  was  introduced  only  about  five  years  back.  Prior  to 
this  research,  all  approaches  used  only  first  generation  wavelets  without  adequate  attention 
to  the  choice  of  mother  wavelets.  Second  generation  scaling  functions  and  wavelets  differ  from 
their  first  generation  counterparts  in  the  property  that  they  are  not  necessarily  translates  and 
dilates  of  one  function.  This  allows  adaptation  to  more  generic  conditions  (such  as  irregular 
samples,  weights,  surfaces  etc.)  while  preserving  some  powerful  and  important  properties 
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of  FGWs  such  as  time-frequency  localization,  multi-resolution  analysis  and  the  existence 
of  fast  transforms.  This  renders  SGWs  useful  in  many  real-life  problems  and  applications 
where  FGWs  cannot  be  readily  applied.  Since  SGWs  sacrifice  the  properties  of  translation 
and  dilation,  they  cannot  be  constructed  using  Fourier  transform  based  methods  and  instead 
rely  on  the  lifting  scheme,  which  is  an  entirely  spatial  domain  construction  technique. 

The  attribute  ‘ blind ’  applied  to  the  superresolution  problem  denotes  the  lack  of  in¬ 
formation  on  the  function(s)  representing  the  blur  in  the  system.  Therefore,  a  HR  re¬ 
construction  algorithm/method  should  include  image  registration,  denoising  and  deblurring 
components  in  addition  to  the  superresolution  module. 

This  research  presents  new  computationally  efficient  algorithms  for  superresolution 
which  offer  better  performance  than  other  current  methods.  Typically,  the  superresolution 
module  is  followed  by  a  separate  denoising  module.  But  in  the  case  of  the  superresolution 
algorithms  presented  in  this  research,  this  module  is  rendered  unnecessary  since  denoising 
is  achieved  simultaneously  with  superresolution,  and  is  hence  implicit  to  the  superresolution 
module.  The  final  step  involves  a  multi-frame  blind  deconvolution  algorithm  which  results  in 
a  denoised  and  deblurred  high  resolution  image.  This  research  improves  on  the  performance 
of  a  popular  and  effective  blind  deconvolution  algorithms  which  extends  directly  to  the 
multi-frame  case. 


3  Summary  of  the  Most  Important  Results  Obtained 

The  use  of  second  generation  wavelets  (SGWs)  to  attain  superresolution  (from  a  captured 
sequence  of  low  resolution  noisy  and  blurred  frames)  with  noise  filtering  has  been  shown  to 
be  possible  without  any  assumption  on  grid  (sampling  lattice)  structure.  The  mathematical 
framework  developed  [1]  is  based  on  nonseparable  two-dimensional  methods  after  demon¬ 
strating  the  feasibility  in  the  separable  case.  The  procedure  allows  the  incorporation  of  the 
more  general  projective  camera  motion  model  into  the  framework,  instead  of  only  displace¬ 
ment  and  rotational  models.  The  primary  goal  of  achieving  blind  superresolution  with  blur 
support  estimation  has  been  realized.  The  results  obtained  proceed  from  a  special  but  useful 
case  of  fast  super-resolution  under  periodic  boundary  condition,  derivation  of  an  optimal 
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threshold  to  resolve  the  tradeoff  between  conflicting  factors  of  satisfactory  blur  removal  and 
noise  reduction  for  a  visually  pleasing  superresolved  image,  to,  subsequently,  multiframe 
blind  superresolution  and  noise  filtering  with,  importantly,  the  accurate  estimation  of  blur 
support. 

3.1  Optimal  Choice  of  Threshold  in  Second  Generation  Wavelet 
Superresolution  [2] 

Wavelet  coefficient  thresholding  has  been  shown  to  be  effective  in  reducing  spatial  domain 
noise  in  the  second  generation  wavelet  based  super-resolution  algorithm.  At  the  same  time, 
thresholding  of  wavelet  coefficients  increases  the  blurring  in  an  image  due  to  loss  of  high 
frequency  information.  Hence,  the  choice  of  threshold  for  noise  reduction  should  be  made  in 
a  manner  so  as  to  achieve  an  optimal  trade-off  between  the  desirable  noise  reduction  and  the 
undesirable  blurring.  In  this  research,  the  effect  of  the  threshold  level  on  reconstructed  image 
quality  has  been  investigated.  In  the  presence  of  blurring  in  the  input  along  with  noise,  the 
choice  of  threshold  plays  an  important  role  in  determining  the  visual  quality  of  the  generated 
high-resolution  image.  The  unsuitability  of  PSNR  as  a  measure  of  image  quality  is  confirmed 
and  a  new  measure  based  on  the  singular  values  of  the  image  is  employed.  Optimal  choice  of 
the  threshold  for  noise  filtering  has  been  coupled  with  adaptive  neighborhoods  for  prediction, 
both  of  which  improve  the  performance  of  the  second  generation  wavelet  superresolution 
(SGWSR)  algorithms. 

3.2  Iterative  PSF  Support  Estimation  for  the  Biggs- Andrews  It¬ 
erative  Blind  Deconvolution  algorithm  [3] 

The  issue  of  blur  in  blind  multiframe  superresolution  has  been  handled,  especially  in  con¬ 
junction  with  an  accurate  estimation  of  unknown  blur  support.  The  iterative  multiframe 
blind  deconvolution  algorithm  considered  by  Biggs- Andrews  requires  exact  knowledge  of  the 
support  of  the  PSF  for  optimal  performance  and  this  is  a  serious  drawback.  Thus,  in  a 
sense,  it  is  not  totally  blind.  A  significant  outcome  of  this  research  is  the  modification  of 
the  Biggs-Andrews  algorithm  so  that  it  incorporates  the  important  step  of  iterative  point 
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spread  function  (PSF)  support  estimation. 


3.3  Moving  Least  Squares  in  Image  Sequence  Superresolution  [4] 

The  Moving  Least  Squares  (MLS)  method  is  a  technique  that  is  likely  to  have  a  place  in 
the  task  of  image  sequence  superresolution  because  of  its  ability  to  interpolate  over  a  grid  of 
irregularly  spaced  data  points  created  by  the  undersampled,  low-resolution,  possibly  blurred 
and  noisy  images  images.  The  basis  functions  chosen  in  approximation  to  get  the  interpolant 
are  usually  multivariate  polynomials,  and  the  highest  total  degree  of  the  polynomial  used  is 
the  order  of  approximation.  A  weighted  least-squares  estimation  procedure  is  used  to  find 
the  coefficients  of  the  basis  function.  The  weight  function  used  is  the  Gaussian  function. 
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ABSTRACT:  The  use  of  second-generation  wavelets  to  attain  super¬ 
resolution  with  noise  filtering  is  described  for  a  captured  sequence  of 
low-resolution  frames  without  any  assumptions  on  grid  (sampling 
lattice)  structure.  The  approach  is  based  on  2D  methods.  The  proce¬ 
dure  allows  the  incorporation  of  the  more  general  projective  camera 
motion  model  into  the  framework,  instead  of  only  displacement  and 
rotational  models.  Several  simulations  that  compare  the  implementa¬ 
tions  of  the  algorithm  presented  here  with  other  related  approaches 
help  illustrate  the  suitability  of  SGWs  (coupled  with  hard  or  soft 
thresholding)  in  the  task  of  image  sequence  superresolution  with 
simultaneous  noise  filtering.  ©  2004  Wiley  Periodicals,  Inc.  Int  J  Imaging 
Syst  Technol,  14,  84-89,  2004;  Published  online  in  Wiley  InterScience  (www. 
interscience.wiley.com).  DOI  10.1002/ima.2001 1 

Key  words:  image  sequence  super-resolution;  second-generation 
wavelets;  high-resolution  reconstruction 


I.  MOTIVATION  AND  INTRODUCTION 

In  this  section  the  suitability  of  second-generation  wavelets  (SGWs) 
for  image  sequence  super-resolution  is  justified  to  be  the  motivating 
factor  behind  the  research  described  here.  This  is  followed  by 
background  material,  with  primary  emphasis  on  the  conditions  im¬ 
portant  in  the  attainment  of  super-resolution,  like  subpixel  displace¬ 
ments  between  adjacent  frames,  because  this  requirement  is  also 
crucial  to  the  procedure  developed  here. 

A.  Motivation.  First-generation  wavelets  (FGWs)  that  are  dilates 
and  translates  of  a  chosen  mother  wavelet  imply  a  regular  sampling 
of  the  data,  are  defined  on  infinite  domains,  and  have  difficulties  in 
dealing  with  boundaries.  At  the  boundary,  where  one  runs  out  of 
data  as  in  the  case  of  a  digital  image  frame  with  a  finite  number  of 
pixels,  zero  padding,  cyclic  and  acyclic  extensions  of  data  have  been 
used  in  the  past.  Special  basis  functions  (translates  and  dilates  of  a 
vector  of  scaling  functions)  (Windmolders  et  al.,  2003)  near  the 
boundary  such  that  all  basis  functions  are  defined  over  the  finite 
extent  array  are  possible  to  use,  but  still  adaptations  near  the  bound¬ 
ary  will  be  required.  Furthermore,  image  sequence  super-resolution 
(high  resolution)  can  be  viewed  in  the  context  of  the  conversion  of 
the  high-resolution  nonuniformly  sampled  raster  (irregular  grid  or 
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nonuniform  sampling  lattice)  generated  from  an  acquired  sequence 
of  low-resolution  frames  to  the  desired  uniformly  sampled  high- 
resolution  grid.  Second-generation  wavelets  (SGWs)  can  not  only 
deal  with  bounded  domains  and  arbitrary  boundary  conditions,  but 
also  irregular  sampling  intervals,  which  are  at  the  heart  of  a  se¬ 
quence  of  low-resolution  images  from  which  image  super-resolution 
is  desired.  FGWs  have  been  used  in  image  sequence  super-resolu¬ 
tion  during  the  last  three  or  four  years  (Bose  and  Lertrattanapanich, 
2004;  Nguyen  and  Milanfar,  2000). 

SGWs  replace  dilations  and  translations  with  an  entirely  spatial 
domain  lifting  scheme  (Sweldens,  1998)  based  on  the  operations  of 
splitting,  prediction,  and  updating.  They  are  suitable  for  use  in 
multidimensions.  In  the  lifting  scheme,  the  forward  transform  is 
implementable  with  low  storage  and  computational  time  costs;  the 
inverse  transform  is  implemented  very  simply  by  essentially  revers¬ 
ing  (updating,  prediction,  and  merging  in  place  of  splitting)  the  order 
of  operations  in  the  forward  transform.  Importantly,  the  lifting 
scheme  adjusts  well  to  boundary  conditions  and  irregular  sampling. 
SGWs  are  a  generalization  of  biorthogonal  wavelets  (traditional 
biorthogonal  wavelets  have  difficulty  dealing  with  boundaries)  and 
form  a  Riesz  basis  for  some  function  spaces  (Sweldens,  1998; 
Sweldens  and  Schroder,  1996;  Vanraes  et  al.,  2002).  SGWs  have 
recently  been  used  to  solve  a  system  of  nonlinear  partial  differential 
equations  with  complicated  boundary  conditions  by  using  a  process 
of  grid  adaptation  and  thresholding  to  control  the  error  in  approxi¬ 
mation  (Vasilyev  and  Bowman,  2000). 

The  area  where  filter  banks  and  wavelets  have  had  the  most 
visible  impact  is  compression.  At  the  other  end  of  the  application 
spectrum,  since  2000,  FGWs  have  been  used  to  attain  super-resolu¬ 
tion  from  image  and  video  sequences.  The  importance  of  the  choice 
of  the  mother  wavelet  and  mother  scaling  function  in  different 
applications  has  been  pointed  out  recently.  For  example,  the  opti¬ 
mum  wavelet  providing  the  maximum  compression  ratio  and  least 
encoding  delay  (compression-speed  tradeoff)  for  audio  signals  from 
stringed  instruments  like  the  sitar  was  found  to  be  the  biorthogonal 
spline  wavelet  whereas  for  audio  signals  from  castanets,  the  Dau- 
bechies  wavelet  was  optimum  (Sathidevi  and  Venkataramani,  2002). 

B.  Introduction.  Techniques  for  enhancing  the  resolution  by  pro¬ 
cessing  the  acquired  images  shifted  by  fraction  of  a  pixel  (subpixel) 
with  respect  to  each  other  are  well  documented.  To  cite  a  few. 
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Poletto  and  Nicolosi  (1999)  developed  a  high-resolution  reconstruc¬ 
tion  algorithm  insensitive  to  high  uncertainties  in  the  relative  sub¬ 
pixel  displacements  of  the  low-resolution  frames  and  to  a  low 
signal-to-noise  ratio.  Shift  information  of  the  sampling  lattices  in 
high-resolution  reconstruction  (image  sequence  super-resolution)  is 
extracted  from  a  given  image  sequence,  static  or  dynamic,  through 
subpixel  accuracy.  In  dynamic  image  reconstruction,  motion  com¬ 
pensation  is  needed  before  subpixel  shift  information  can  be  ex¬ 
tracted  from  the  sequence.  A  survey  of  several  subpixel  accuracy 
image  registration  algorithms  for,  especially,  machine  vision  algo¬ 
rithms  is  available  in  Kim  and  Su  (1993).  In  Su  and  Kim  (1994),  for 
each  block  image  in  a  reference  frame,  motion  estimation  and 
subpixel  registration  are  performed  against  adjacent  frames,  prior  to 
high-resolution  reconstruction  from  dynamic  image  sequences. 

Algorithms  have  been  proposed  that  use  subpixel  shifts  between 
a  sequence  of  captured  images  [low-resolution  (LR)  frames]  to 
estimate  a  single  high-resolution  (HR)  image,  called  a  super-re- 
solved  image.  This  process  has  been  referred  to  as  microscanning 
(Gillett,  et  al.,  1995).  Such  subpixel  shifts  may  occur  due  to  random 
variation  between  the  objects  in  a  scene  and  camera  motion.  The 
sampling  of  the  scene  is  thus  nonuniform,  and  the  shifts  are  com¬ 
bined  so  that  an  image  that  is  effectively  sampled  at  a  higher  rate 
than  the  individual  frames  is  constructed.  Improved  resolution  from 
subpixel  shifted  pictures  was  also  obtained  (Ur  and  Gross,  1992)  by 
using  a  multichannel  generalized  sampling  theorem.  The  need  for 
subpixel  displacements  was  apparent  in  the  recursive  algorithms 
developed  for  HR  reconstruction  from  LR  frames  in  the  presence  not 
only  of  observation  error  but  also  error  in  estimation  of  the  unknown 
shift  parameters  (Bose  et  al.,  1993).  For  subpixel  motion  estimation 
the  work  of  Schultz  et  al.  (1998)  is  also  relevant.  Bose  and  Boo 
(1998)  presented  a  mathematical  model  for  constructing  an  HR 
image  from  LR  images  acquired  through  a  multisensor  array.  Sub¬ 
pixel  shifts  were  needed  in  this  model  and,  therefore,  in  subsequent 
research  that  used  the  model,  (e.g.,  Chan  et  al.,  2003;  Ng  et  al., 
2002). 

In  the  early  work  on  wavelet  superresolution,  the  important 
problem  of  mother  wavelet  selection  was  neglected  (Nguyen  and 
Milanfar,  2000).  The  finitely  supported  families  evaluated  for  the 
task  of  image  sequence  super-resolution  [discussed  in  Lertrattan- 
apich  (2003)]  for  optimum  selection  of  mother  wavelet  included 
Haar,  Daubechies,  Symlets,  Coiflets,  and  B-splines.  The  common 
desirable  properties  of  mother  scaling  function  and  wavelet — in¬ 
cluding  finite  small  support  (to  reduce  computational  chore),  explicit 
and  simple  expression,  symmetry  (or  linear  phase),  orthogonality  or 
biorthogonality,  high  regularity  (smoothness),  low  approximation 
error,  and  high  time — frequency  localization — are  best  satisfied  by 
B-splines  among  the  classes  considered. 

II.  OBJECTIVE 

The  basic  idea  of  second-generation  wavelets  (SGWs)  is  to  sacrifice 
translation  and  dilation  in  order  to  construct  wavelets  adapted  to 
general  settings,  but  to  simultaneously  preserve  many  powerful 
properties  of  first-generation  wavelets  (FGWs)  such  as  time — fre¬ 
quency  localization,  orthogonality  or  biorthogonality,  and  fast  trans¬ 
forms.  It  has  been  noted  quite  recently  that  SGWs  are  inherently 
more  suited  for  image  super-resolution.  In  Bose  et  al.  (2004),  the  HR 
reconstruction  with  noise  reduction  was  restricted  to  semi-regular 
grids  (sampling  lattices)  [tensor  product  of  ID  irregular  grids  (sam¬ 
pling  lattices)];  consequently,  the  LR  frames  were  displaced  only  by 
translations  from  a  reference  frame,  as  approximated  in  far-held 
imaging — e.g.,  satellite  imaging  or  aerial  photography. 


In  the  research  reported  here,  attention  is  focused  on  achieving 
super-resolution  from  registered  noisy  LR  frames  without  any  as¬ 
sumptions  on  grid  (sampling  lattice)  structure:  that  is,  the  procedure 
will  be  able  to  handle  arbitrary  irregular  grids  (sampling  lattices), 
although  the  subpixel  displacement  assumption,  which  is  justified  in 
super-resolution  algorithms,  will  be  maintained.  A  consequence  of 
the  use  of  arbitrary  grids  will  be  the  ability  to  incorporate  the  more 
general  projective  camera  motion  model  into  the  framework.  The 
ability  to  handle  arbitrary  sampling  lattices  results  from  the  use  of 
2D  prediction  and  update  operators  as  opposed  to  the  tensor  prod¬ 
uct  of  ID  operators  used  by  Bose  et  al.  (2004).  Furthermore,  in  the 
2D  case,  the  definition  of  odd  and  even  samples  in  ID  cannot  be 
directly  extended.  This  subsequently  results  in  the  inability  to  define 
levels  of  a  multiresolution  decomposition  in  the  dyadic  sense.  Thus, 
in  the  SGW  setting,  only  one  coefficient  per  scale  is  chosen  for 
prediction  (Vanraes  et  al.,  2002).  The  use  of  genuine  2D  methods 
gives  better  results,  as  seen  from  simulation  results  presented  below. 

III.  PREPROCESSING  PROCESS  FOR  SUPER¬ 
RESOLUTION  ALGORITHM 

A.  Singularity  Prevention.  In  the  super-resolution  problem,  af¬ 
ter  registering  all  samples  in  LR  images  into  a  HR  grid  (sampling 
lattice),  it  is  possible  that  two  samples  or  more  are  either  trans¬ 
formed  to  the  same  location  or  assumed  to  be  the  same  because  of 
machine  imprecision.  In  some  super-resolution  algorithms,  certain 
irregular  sampling  structures  in  a  HR  grid  (sampling  lattice)  that  can 
lead  to  a  singular  problem  were  identified  by  Kim  and  Bose  (1990). 
Therefore,  here  also,  the  LR  frames  that  give  rise  to  such  singular 
sampling  structures  are  detected  and  eliminated  before  the  super¬ 
resolution  algorithm  is  implemented. 

B.  Bad  Frame  Elimination.  Another  possible  problem  arises 
because  some  LR  frames  in  the  sequence  may  be  too  badly  cor¬ 
rupted  by  additive  noise  or  severely  blurred  (or  both)  to  be  useful  in 
processing  when  compared  with  other  frames  in  the  sequence.  Such 
frames  are  called  “bad  frames.”  These  bad  frames  should  be  dis¬ 
carded  from  the  sequence  because  they  usually  deteriorate  the  re¬ 
sulting  HR  image.  A  rationale  for  identifying  the  bad  frames  has 
been  advanced  recently  by  Lertrattanapanich  (2003).  Because  most 
parts  of  LR  frames  usually  contain  the  same  region  of  interest  (ROI), 
the  variances  of  these  LR  frames  in  the  sequence  should  be  more  or 
less  the  same  unless  they  are  bad  frames.  Note  that  the  sample 
variance  of  of  the  Mi  M  X  N  LR  frame  fk[i,  j\  used  here  is  defined 
as 

1  M-l  N- 1 

°i  =  MN  _  t  X  S  (/*[*>  j\  ~  m*)2.  (i) 

i=0  j= 0 

where 

1  M-l  N- 1 

lb  —  ^  2  fk\_f  j \  (2) 

i=0  j= 0 

is  the  sample  mean  of  the  Mi  LR  frame  fk[i,  /]. 

Intuitively,  the  bad  frames  that  are  badly  corrupted  by  additive 
noise  have  higher  variance  than  acceptable  LR  frames.  In  contrast, 
the  variances  of  the  bad  frames  that  are  severely  blurred  by  some 
PSFs  are  lower  than  those  of  acceptable  LR  frames.  This  expectation 
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has  been  substantiated  by  numerous  simulations  (Lertrattanapanich, 
2003)  and  has  been  useful  in  the  identification  and  subsequent 
removal  of  bad  frames. 

IV.  SUPER-RESOLUTION  ALGORITHM  USING  SGWs 

In  this  section,  an  algorithmic  description  of  the  process  for  obtain¬ 
ing  a  super-resolved  image  from  a  set  of  LR  images  based  on  SGWs 
and  the  lifting  technique  (Sweldens,  1998)  is  presented.  Lifting  is  a 
general  technique  that  can  be  used  for  the  construction  of  FGWs  and 
SGWs.  FGWs  can  also  be  constructed  using  filter  banks,  but  SGWs 
can  only  be  constructed  using  the  lifting  technique.  The  algorithm 
presented  represents  a  new  and  seemingly  natural  approach  to  super- 
resolution.  In  addition,  noise  filtering  is  simultaneously  imple¬ 
mented. 

Assumptions:  The  main  assumption  made  is  that  the  LR  frames 
have  subpixel  displacements  with  respect  to  a  chosen  reference  LR 
frame,  which  is  a  valid  assumption  in  the  super-resolution  frame¬ 
work  (see  e.g.,  Schultz  et  al.,  1998;  Ur  and  Gross,  1992) 

A  secondary  assumption  is  that  the  LR  frames  are  either  already 
registered  or  their  registration  matrices  (projective  camera  motion 
model  parameter  matrices)  are  known  by  the  use  of  standard  pro¬ 
cedures  (Lertrattanapanich,  2003;  Lertrattanapanich  and  Bose,  2002; 
Mann  and  Picard,  1997). 

Let  LRfm,  n],  i  =  1,  2,  ....  r,  ....  K  represent  the  LR  frames, 
where  LRr  denotes  the  reference  LR  frame.  The  coordinate  system 
considered  is  with  respect  to  LRr.  Let  J\  m,  ri\  denote  the  registered 
image  obtained  using  a  preselected  subset  Kl  C  { 1,  2,  .  .  .  ,  K)  of  the 
available  noisy  frames.  Let  s[m,  n]  denote  the  super-resolved  image 
that  exists  on  a  regular  HR  grid,  the  coordinates  of  which  are  defined 
by  the  rows  of  a  matrix  H.  The  order,  N,  of  subdivision  is  determined 
by  the  bivariate  polynomial  of  partial  degree  (N  —  1)  in  each 
variable  used  for  prediction  and  can  be  shown  (Sweldens,  1998)  to 
correspond  to  the  number  of  vanishing  moments  of  the  dual  wavelet 
(counterpart  of  scaling  function  in  FGW),  <p.  The  number  of  van¬ 
ishing  moments  of  the  primal  wavelet,  if/,  N,  is  chosen  in  the  update 
step.  A  closed  form  expression  cannot  be  given  for  the  dual  and 
primal  wavelets  (unlike  in  some  FGWs  like  B-spline  based)  because, 
in  the  general  second-generation  setting,  wavelets  are  custom  built 
to  fit  the  given  samples,  and  hence  are  arbitrary.  In  the  simulation 
results  presented  here,  N  =  2,  N  =  1,  because  a  1-ring  neighborhood 
is  used  for  prediction  (note  that  for  higher  values  of  N  and  N,  a  larger 
neighborhood  will  be  required  to  ensure  a  full  rank  overdetermined 
system,  used  here,  in  the  prediction  operation  described  subse¬ 
quently).  Reasons  for  the  choice  of  a  1-ring  neighborhood  in  place 
of  a  1-ring  neighborhood  with  flaps  are  presented  in  the  algorithm 
given  below.  Other  higher-order  neighborhoods,  such  as  the  2-ring 
neighborhood,  are  supersets  of  the  1-ring  neighborhood  with  flaps, 
and  hence  cannot  be  used  for  the  same  reasons  as  given  for  the 
1-ring  neighborhood  with  flaps. 

The  least  squares  prediction  at  a  generic  point  [m0,  n0\  can  be 
viewed  as  the  value  at  [m0,  n0 ]  of  the  surface  that  fits,  in  a  least 
squares  sense,  the  points  in  N  ,  the  1-ring  2D  neighborhood  of 
[m0,  n0].  Let  in  =  (mh  .  .  . ,  mx]T  and  n  =  [nlt  .  .  .  ,  represent 
the  vectors  associated  with  the  points  belonging  to  N  ,  where  JV 
is  the  number  of  points  in  N  .  Let  s  represent  the  corresponding 
vector  of  scaling  function  coefficients.  The  model  for  the  surface  (a 
plane  in  this  case  since  N  =  2)  is  given  by 

s  =  a0  +  apn  +  a2n  +  e, 

where  e  is  the  error  in  the  fit  with  mean  E{e)  =  0.  The  minimum 
norm  least  squares  estimate  of  a  =  [a0  oq  a2]  is  given  by 


a  =  (PTPr‘PTs, 

where  the  matrix  P  =  [1  m  n]  and  1  is  a  vector  of  l’s.  The  prediction 
at  a  generic  point  [m0,  «0]  is  then  computed  as  the  dot  product 

s0  =  [1  mo  «o]  •  a. 

The  detail  coefficient  that  replaces  the  original  value  s0  at  [m0,  n0] 
is  then  computed  as 


do  —  So  So-  (3) 

The  update  operator  can  be  described  by  the  vector  u  = 
[iq,  .  .  .  ,  uM]  where  each  of  the  u-  s  is  computed  so  that  its  €2-norm 
is  minimized  using  the  equation  (Delouille  et  al.,  2003;  lansen  et  al., 
2001), 


Ui  = 


AA*0 
(A*)2  ’ 


(4) 


where  Aq  is  the  support  area  of  (f>  around  the  point  [w0,  n0\  and  A,- 
denotes  the  corresponding  area  around  point  i  £  N  after  the  grid 
has  been  coarsened  by  the  removal  of  the  predicted  point  [w0,  n0]. 
Subsequently,  each  scaling  function  coefficient,  sh  is  updated  ac¬ 
cording  to  a  set  of  equations  associated  with  all  the  detail  coeffi¬ 
cients  predicted  using  s t  A  typical  member  of  this  set  associated 
with  the  generic  point  [m0,  n0]  is 


x 

do  =  s0  -  X  skuk.  (5) 

i-=i 


A.  Algorithm. 

•  If  Tim,  n]  is  not  given,  obtain  f[m,  n\  using  the  known  regis¬ 
tration  matrices  of  LRfm,  n],  i  £  on  the  coordinate  system 
IR  X  [R  defined  with  respect  to  LR,.  The  grid  (sampling  lattice) 
(on  which  f[m,  n]  exists)  so  obtained  is  an  irregularly  sampled 
grid  (sampling  lattice). 

•  Define  the  required  high-resolution  grid  (sampling  lattice)  in 
the  coordinate  system  defined  above,  that  is,  determine  the 
matrix  H.  Split  H  into  two  sets  of  points,  one  set,  Hu,  com¬ 
prised  of  points  whose  values  are  to  be  updated  and  retained  in 
the  subsequent  steps  of  the  algorithm,  and  another  set,  Hp, 
which  is  to  be  predicted  and  eliminated  from  the  coarse  rep¬ 
resentation  offlm,  n\  developed  below. 

•  Compute  the  Delaunay  triangulation  (Okabe  et  al.,  1992)  of  the 
points  defining /[m,  n\. 

•  For  each  point  in  the  set  Hp,  say  [m0,  n0]: 

-Determine  the  points  in  the  1-ring  neighborhood  (Figure  1), 
Nm  ,  to  the  extent  possible.  The  use  of  a  1-ring  neighbor¬ 
hood  with  flaps  for  prediction  was  documented  by  Delouille 
et  al.  (2003).  But  this  was  not  employed  because  the  assump¬ 
tion  of  subpixel  displacement  of  the  LR  frames  tends  to  result 
in  a  triangulation  with  flap  points  that  are  usually  nearly 
collinear  with  the  points  in  the  1-ring  neighborhood.  This 
may  lead  to  ill  conditioning  of  the  prediction  matrix  which, 
consequently,  causes  instability  in  the  implementation  of  the 
prediction  operator. 
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1  -ring  neighborhood 


-Compute  the  prediction  at  [ m0 ,  «0]  by  means  of  a  least 
squares  predictor,  described  above  [boundary  points  are  han¬ 
dled  differently  as  in,  for  example],  Delouille  et  al.  (2003), 
using  points  in  N  .  Subsequently,  replace  the  scaling  func¬ 
tion  coefficient  at  [m0,  n0]  by  the  detail  coefficient  in  Eq.  (3). 
-Determine  A^,  the  support  area  of  the  scaling  function  around 
[m0,  n()  | .  This  can  be  determined  as  the  sum  of  one-third  the 
area  of  the  triangles  with  common  vertex  [ m0 ,  n0], 

-Delete  [m0,  n0]  from  the  triangulation  resulting  in  a  coarser 
grid  (sampling  lattice) 

-Compute  the  corresponding  areas,  Ah  k  =  1,  .  .  . ,  X,  of  the 
cells  around  points  belonging  to  N  in  the  coarsened  grid. 
-Compute  the  minimum  norm  update  operator  as  in  Eq.  (4) 
and  update  the  scaling  function  coefficients  at  points  belong¬ 
ing  to  Nmono,  as  in  Eq.  (5). 

A  second-generation  2D  wavelet  surface  has  now  been  defined 
on  the  irregularly  sampled  registered  grid  (sampling  lattice). 
The  detail  coefficient  values  at  the  points  on  the  HR  grid 
(sampling  lattice)  can  now  be  computed  either  by  ignoring  the 
irregularity  of  the  grid  and  assigning  the  value  of  the  closest 
irregular  point  (no  significant  loss  of  accuracy  due  to  subpixel 
displacement  assumption)  or  by  resampling  the  wavelet  (not 
scaling  function). 


•  To  achieve  noise  reduction,  the  wavelet/detail  coefficients 
computed  are  subjected  to  thresholding  (Jansen,  2001;  Vanraes 
et  al.,  2002). 

•  The  inverse  lifting  procedure  is  now  applied  on  the  HR  grid 
(sampling  lattice)  to  obtain  the  super-resolved  image  s\m,  n\. 

V.  SIMULATION  RESULTS 

Results  obtained  by  the  use  of  the  SGW-based  algorithm  described 
in  the  preceding  section  are  presented  and  discussed  below.  The  LR 
frames  used  in  the  simulation  were  of  size  40  X  40  pixels  and  were 
generated  using  an  80  X  80  image  by  downsampling  and  low-pass 
filtering.  Further,  the  LR  frames  are  all  displaced  with  respect  to 
each  other,  but  unlike  the  results  presented  by  Bose  et  al.  (2004)  (in 
which  the  displacements  were  restricted  to  only  translations),  the 
displacements  follow  the  more  general  projective  model.  The  LR 
frames  were  also  corrupted  by  AWGN  of  variance  0.1.  As  men¬ 
tioned  earlier,  in  the  simulations  presented  in  this  section,  the 
number  of  vanishing  moments  of  the  dual  and  primal  wavelets,  N 
and  N  respectively,  are  2  and  1. 

Figure  2  shows  the  results  from  one  simulation.  In  this  case,  six 
LR  frames  were  used  and  a  resolution  improvement  factor  of  2  in 
each  dimension  was  chosen.  It  is  noted  that  the  use  of  higher  number 
of  LR  frames  leads  not  only  to  higher  resolution  factors,  but  also  to 
improved  noise  reduction.  A  sample  LR  frame  is  shown  in  Figure 
2(a)  and  the  original  image  is  shown  in  Figure  2(b).  The  SGW-based 
high-resolution  reconstructions  using  hard  and  soft  thresholding  are 
shown  in  Figures  2(c)  and  2(d),  respectively.  For  comparison,  four 
other  high-resolution  reconstructions  based  on  surface  approxima¬ 
tion  are  given  in  Figures  2(e) — (h).  Note  that  direct  interpolation 
(bilinear  and  bicubic)  of  a  single  LR  frame  gives  the  visually  worst 
results,  although  their  PSNR  value  is  comparable  to  the  others.  The 
Delaunay  triangulation-based  methods  give  much  better  results  but 
cannot  simultaneously  achieve  noise  filtering. 

Results  for  another  image  are  shown  in  Figure  3.  As  in  the 
previous  result,  six  LR  frames  are  used  and  the  resolution  improve¬ 
ment  factor  is  set  to  2  in  each  dimension.  Panel  (a)  shows  a  sample 
LR  frame.  The  original  image  is  shown  in  Figure  3(b).  The  SGW- 


SGWSR  (hard 
thresholding) 
21.6336  dB 


SGWSR  (soft 
thresholding) 
22.4785  dB 


bilinear  interpolation 
of  LR  frame  2 
21.0683  dB 


bilinear  interpolation  bicubic  interpolation 
after  DT  of  LR  frame  2 
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bicubic  interpolation 
after  DT 
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Figure  2.  Simulation  result  1  with  PSNR  in  dB. 
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Sample  LR  frame 


Original  image 


SGWSR  (hard 
thresholding) 
23.1867  dB 


SGWSR  (soft 
thresholding) 
22.8185  dB 


bilinear  interpolation 
of  LR  frame  2 
21.7991  dB 


bilinear  interpolation 
after  DT 
22.1385  dB 


bicubic  interpolation 
of  LR  frame  2 
21.5806  dB 


bicubic  interpolation 
after  DT 
21.8509  dB 


Figure  3.  Simulation  result  2  with  PSNR  in  dB. 


based  HR  reconstructions  with  noise  reduction  achieved  using  hard 
and  soft  thresholding  of  wavelet  coefficients  are  shown  in  Figures 
3(c)  and  3(d),  respectively.  Similar  to  the  previous  result  presented, 
surface  interpolation  based  reconstructions  are  given  in  Figures 
3(e)— (h). 

In  both  the  results  presented,  a  point  to  be  noted  is  the  retention 
of  a  reasonably  high  degree  of  sharpness  in  the  noise-filtered  SGW- 
based  reconstructions.  This  is  in  contrast  to  other  noise  filtering 
techniques  such  as  averaging  and  low-pass  filtering,  which  typically 
result  in  smoothing  of  edges. 

To  enable  comparison  and  differentiation  of  the  techniques  pre¬ 
sented  in  this  article  and  in  Bose  et  al.  (2004),  results  generated  by 
both  the  methods  applied  to  the  same  data  set  are  presented  in  Figure 


4.  In  this  simulation,  the  displacements  of  the  LR  frames  were 
restricted  to  translations  [required  to  run  the  algorithm  given  by 
Bose  et  al.  (2004)]  and  four  LR  frames  were  used.  The  LR  frames 
were  corrupted  by  AWGN  of  a  very  low  variance  of  0.05.  In  both 
methods,  N  =  2,  N  =  1.  Panels  (a)  and  (b)  give  a  sample  LR  frame 
and  the  original  image  respectively.  Figures  4(c)  and  4(d)  show  the 
reconstructions  obtained  by  treating  the  registered  image  j\m,  n  |  as 
the  tensor  product  of  two  ID  vectors.  On  treatment  of  f[m,  n  |  in  the 
2D  sense,  the  HR  images  obtained  are  given  in  Figures  4(g)  and 
4(h).  It  is  important  to  observe  that  although  the  two  sets  of  results 
have  very  close  PSNR  values,  the  results  obtained  using  the  algo¬ 
rithm  presented  here  are  much  sharper.  For  comparison,  the  results 
obtained  using  bilinear  and  bicubic  surface  interpolation  after 


1  -D  SGWSR  (hard 
thresholding) 
28.0125  dB 


1  -D  SGWSR  (soft 
thresholding) 
28.1108  dB 


bilinear  interpolation 
after  DT 
26.6031  dB 


bicubic  interpolation 
after  DT 
26.1083  dB 


2-D  SGWSR  (hard 
thresholding) 
28.0055  dB 


2-D  SGWSR  (soft 
thresholding) 
28.0260  dB 


Figure  4.  Comparison  of  ID  (Bose  et  al.,  2004)  and  2D  SGWSR  algorithms  with  PSNR  in  dB. 
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Delaunay  triangulation  are  given  in  Figures  4(e)  and  4(f),  respec¬ 
tively.  Clearly,  both  the  SGWSR  methods  outperform  the  surface 
interpolation  methods. 

VI.  CONCLUSIONS 

A  framework  for  achieving  image  sequence  super-resolution  simul¬ 
taneously  with  noise  filtering  has  been  developed  using  SGWs 
coupled  with  hard  or  soft  thresholding.  The  procedure  works  well 
when  the  noisy  images  acquired  by  sensors  are  not  severely  blurred. 
A  topic  of  future  research  is  the  embedding  of  a  technical  device  into 
the  procedure  to  adequately  overcome  the  degrading  effects  of  blur, 
when  present  to  a  significant  extent.  The  problem  is  reminiscent  to 
that  faced  in  earlier  works  on  super-resolution.  For  example,  in  the 
early  DFT-based  procedure  (Kim  and  Bose,  1990),  interpolation  and 
noise  filtering  were  simultaneously  implemented  and  it  was  later 
shown  that  blur  distortions  in  the  input  images  can  be  compensated 
with  regularization  during  reconstruction  (Kim  and  Su,  1993).  In 
recent  super-resolution  work,  multiframe  blur  identification  fol¬ 
lowed  by  deblurring  has  been  done  in  a  separate  module  from  the 
super-resolution  and  noise  filtering  modules  (Lertrattanapanich  and 
Bose,  2002). 

The  procedure  depends  on  subpixel  accuracy  image  registration. 
For  dynamic  image  sequences,  motion  compensation  is  necessary 
prior  to  that  and,  although  the  algorithm  presented  here  applies  to 
arbitrary  nonuniform  shift  grids  (sampling  lattices),  its  robustness  to 
errors  in  estimation  of  motion  parameters  is  another  topic  of  future 
research. 

The  main  advantage  of  the  procedure  is  that  the  speed  of  imple¬ 
mentation  provided  by  the  lifting  technique  (surface  approximation 
methods  (Lertrattanapanich  and  Bose,  2002)  may  be  faster  but  do 
not  incorporate  noise  filtering),  the  adaptation  to  arbitrary  nonuni¬ 
form  sampling  lattice  (especially  relevant  in  dynamic  image  se¬ 
quences),  the  absence  of  a  priori  assumptions  on  boundary  condi¬ 
tions  (zero,  periodic,  Neumann  have  given  progressively  better 
results  over  the  years  but  make  assumptions  that  are  unnecessary 
here),  and  the  independence  from  proper  choice  of  mother  wavelets 
and  scaling  functions  which  makes  SGW-based  super-resolution 
methods  potentially  more  suitable  in  multimedia  applications. 
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Abstract — Wavelet  coefficient  thresholding  is  effective  in  re¬ 
ducing  spatial  domain  noise  in  wavelet-based  super-resolution 
algorithms.  Here,  the  effect  of  the  threshold  level  on  reconstructed 
image  quality  in  second-generation  wavelet  super-resolution  is 
investigated.  The  choice  of  optimal  threshold  involves  a  trade-off 
between  noise  filtering  and  blurring  introduced  by  thresholding. 
A  measure  based  on  the  singular  values  of  the  image  matrix  is 
employed  as  a  reliable  gauge  of  generated  high-resolution  image 
quality. 

Index  Terms — Noise  filtering,  second-generation  wavelets, 
super-resolution,  thresholding. 

I.  Introduction 

In  this  letter,  the  term  “super-resolution”  refers  to  algorithms 
that  produce  a  high-resolution  (HR)  image  by  combining 
information  from  a  captured  sequence  of  low-resolution  (LR) 
frames.  Due  to  various  factors  like  imperfections  in  the  ac¬ 
quisition  device,  limited  resolution  of  the  physical  sensing 
elements,  motion,  and  medium  turbulence,  the  LR  frames 
are  blurred  and  noisy.  Hence,  super-resolution  algorithms 
commonly  include  noise  filtering  and  deblurring  modules. 
The  focus  of  this  letter  is  to  improve  the  performance  of 
the  second-generation  wavelet  super-resolution  (SGWSR)  al¬ 
gorithms  given  initially  in  [1]  and  in  a  more  general  tow¬ 
dimensional  (2-D)  setting  in  [2],  though  some  of  the  results 
presented  apply  to  wavelet  denoising  in  general.  The  aim  is  to 
remove  as  much  of  the  corrupting  noise  as  possible  without 
adversely  affecting  the  reconstructed  image  quality  due  to  blur 
introduced  in  the  process.  Though  there  are  many  techniques 
for  denoising  images  (for  example  [3],  [4],  [5]  and  [6]), 
some  of  which  consider  optimal  thresholding,  the  denoising 
technique  based  on  wavelet  coefficient  thresholding  [7]  is 
selected  here  because  this  allows  simultaneous  noise  filtering 
and  superresolution  thereby  reducing  computation  cost  in  the 
overall  framework.  Most  of  the  earlier  methods  dealing  with 
optimal  thresholding  of  wavelet  coefficients  adopt  a  quality 
metric  based  on  mean-squared  error  (MSE),  which  does  not 
necessarily  yield  the  best  visual  quality.  The  optimal  threshold 
for  denoising  based  on  wavelet  coefficient  thresholding  is 
obtained  in  this  paper  in  terms  of  the  visual  quality  defined  by 
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a  metric  that  matches  the  human  visual  system  (HVS)  better 
than  other  metrics  [8].  The  SGWSR  methods  simultaneously 
incorporate  noise  filtering  by  hard/soft  thresholding  of  wavelet 
coefficients  defined,  respectively,  by  the  following  expressions 

[7]: 


dh 


d,  d  >  ( 

0,  d^C  ’ 


d-C,  d>  C 

0,  cKC 


0) 


where  d  represents  the  SGW  coefficients  before  thresholding 
and  (  is  the  threshold.  Soft  thresholding,  known  to  yield  better 
results  [1],  [2],  is  used  in  this  paper.  The  algorithms  in  [1], 
[2]  did  not  investigate  the  optimal  choice  of  threshold,  £.  It  is 
important  to  note  that  thresholding  in  the  SGW  setting  is  much 
more  challenging  than  in  the  first-generation  wavelet  (FGW) 
case  since  small  coefficients  may  carry  important  information 
which  is  essential  for  a  stable  inverse  transform  [9].  The  trade¬ 
off  between  blur  introduced  by  thresholding  and  input  noise 
mandates  the  need  to  find  an  optimal  threshold  as  explained  in 
Section  II.  Section  III  substantiates  why  peak  signal-to-noise 
ratio  (PSNR)  and  the  universal  image  quality  index  Q  [10]  are 
unsuitable  as  measures  of  visual  quality  of  an  image.  An  image 
quality  measure  based  on  the  singular  values  of  the  image 
matrix  is  then  used  to  assess  the  effects  of  thresholding  on 
reconstructed  image  quality  in  Section  IV.  Finally,  conclusions 
are  presented  in  Section  V. 


II.  Thresholding  SGW  Coefficients 

A  set  of  sub-pixel  shifted  frames,  { /fc } ,  k  =  1, 2, . . . ,  K  (of 
size  ni  x  n2)  of  the  same  object/scene,  generated  with  small 
variation  in  camera  motion  and  position,  and  captured  under 
the  same  imaging  conditions  should  have  approximately  the 
same  variance,  <j\,  and  mean,  /i^,  given  by 
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When  zero-mean  additive  white  noise  is  present,  the  variance 
of  the  corrupted  frame,  will  be  given  by 


~2  2,2 
ak  ~  ak  +  V  > 


(4) 


where  lj2  is  the  variance  of  the  noise.  A  severe  additive  noise 
is  one  whose  variance  is  high.  Therefore,  a  frame  which  is 
corrupted  by  a  severe  noise  will  have  a  variance  much  higher 
than  that  of  a  frame  with  a  high  SNR.  Blurring  an  image  with 
a  point  spread  function  (PSF)  h\p,  q ],  of  support  to  1  x  to2,  is 


IEEE  SIGNAL  PROCESSING  LETTERS,  VOL.  12,  NO.  11,  NOVEMBER  2005 


2 


viewed  as  low-pass  filtering  of  the  image.  The  variance  of  the 
blurred  frame,  fk,  can  be  proved  to  be 
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h\p,q]h[r,s]Cfk[(p,r);  (g,s)],  (5) 


where  Cfk  is  the  autocovariance  of  the  original  image,  fk  ■ 
If  the  PSF  is  assumed  to  sum  to  unity  over  its  support,  then 
since  |C/JfeO;(g,s)]|  <  Cfk[(p,p);(q,q)\  =  a2k. 
Consequently,  blur  tends  to  reduce  the  variance  unlike  noise, 
which  as  seen  from  (4),  increases  the  variance. 

The  opposing  effects  of  blur  and  noise  on  the  variance  also 
hold  when  the  SGW  coefficients  are  thresholded.  Similar  to  the 
FGW  setting,  the  noise  corrupting  an  image  is  carried  over  into 
the  SGW  coefficients  on  application  of  the  forward  transform. 
This  is  readily  seen  from  the  expression  for  the  computation 
of  the  SGW  coefficients  (based  on  lifting  [11])  given  by 


d-i[i,j\  =  (s0[*>  j]  +  n[hj])  ~  P{N[i,j]},  (6) 

where  so[i,j]  is  the  pixel  value  at  location  [i,j\,  n [i,j] 
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to  be  indicators  of  better  visual  quality.  However,  it  is  known 
that  PSNR  does  not  match  the  HVS  [12];  that  is,  higher 
PSNR  does  not  always  imply  better  visual  quality.  Studies 
have  shown  that  the  HVS  is  more  sensitive  to  blurring  than 
to  noise.  Another  image  quality  measure  Q  proposed  in  [10], 
though  more  sensitive  to  blur  than  PSNR,  was  shown  by  our 
simulations  to  be  not  ideal  for  determining  £0pt,  the  optimal 
threshold,  consistent  with  expectations  from  results  in  [8],  Q 
is  defined  as  [10] 


Q  = 


(9) 


where  [if  and  [ij  are  the  means,  cr 2  and  a 2  are  the  variances 
of  the  original  and  the  observed  images,  respectively,  and 
is  the  cross-covariance  between  the  original  and  the  observed 
images.  The  plots  of  PSNR  and  Q  versus  the  threshold  for  an 
SGWSR  reconstruction  are  shown  in  Fig.  2.  The  input  to  the 
SGWSR  algorithm  is  a  set  of  six  registered  LR  frames,  fk, 
generated  from  the  original  image,  /,  as 


fk  =  DHTkf  +  77-fc,  k  =  1, . . . ,  6,  (10) 


where  Tk  is  the  projective  transformation,  H  is  the  blur  matrix 
associated  with  a  5  x  5  truncated  Gaussian  PSF  (note  that 
this  blur,  which  affects  the  input  frames,  is  different  from  the 
blur  introduced  due  to  thresholding  in  the  reconstructed  HR 
images),  D  is  the  downsampling  matrix,  and  rik  is  additive 
white  Gaussian  noise  (AWGN)  of  variance  0.01.  Despite  the 


Fig.  1 .  Effect  of  threshold  on  input  noise  and  blur  introduced  by  thresholding 

represents  random  white  noise,  d_i[i,j]  is  the  computed  SGW 
coefficient,  and  P  is  the  prediction  operator  that  operates  on 
the  specified  neighborhood  of  point  A  detailed 

description  of  the  process  of  prediction  for  SGWSR  is  given 
in  [2],  The  prediction  operator,  P,  fails  to  predict  the  white 
noise  which  propagates  to  the  SGW  coefficient  d-\[i,j]  in  (6). 
Thus,  thresholding  of  the  wavelet  coefficients  reduces  the  noise 
in  the  reconstruction.  However,  high  wavenumber  information 
in  the  SGW  coefficients  may  be  eliminated  by  thresholding 
leading  to  possibly  significant  blurring  of  the  reconstructed 
image.  Thus,  a  need  for  an  optimal  choice  of  the  threshold 
exists.  Fig.  1  illustrates  the  above  points. 

III.  Image  Quality  Measures 
PSNR  and  MSE  are  the  most  commonly  used  measures  for 
comparative  image  quality.  They  are  defined  for  normalized 
ni  x  n 2  images  as 

mse  =  (7) 

n\ri2 

PSNR=201o&”bk)-  f8) 

where  /  and  /  represent  the  original  image1  and  the  observa¬ 
tion,  respectively.  Higher  PSNR  and  lower  MSE  are  expected 

[n  real- world  problems,  the  original  image  is  usually  inaccessible,  which 
motivates  Section  IV 


Fig.  2.  Effect  of  threshold  on  PSNR  and  Q  in  SGWSR 

increased  blurring  and  visual  quality  degradation  with  higher 
thresholds  (as  observed  in  Fig.  1),  both  PSNR  and  Q  measures 
are  inadequate  and  not  sufficiently  sensitive  to  this  visual 
quality  change.  In  Fig.  2,  though  the  curve  appears  to  be 
perfectly  flat  before  and  after  the  threshold  marked  as  Copt, 
there  is  a  slight  increase  in  PSNR  and  Q  values  at  Copt  (see 
insets  in  Fig.  2)  after  which  they  remain  constant.  Flattening 
of  the  curve  occurs  since  all  wavelet  coefficients  fall  within 
the  threshold  as  it  is  increased  beyond  Copt-  When  considering 
such  flat  regions,  the  minimum  possible  threshold  is  selected 
US  Copt  - 
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Another  measure  of  image  quality  based  on  the  singular 
values  of  the  image  matrix  was  defined  as  [8] 


Msvd 


E^E  \Dj-Dmed\ 
ki  x  fc2 


Dj  = 


\|B»*  «)’, 


(ID 


where  k\  =  ni/p,  fc2  =  112/ p,  Dmed  is  the  median  of  the  Dj’s 
and  Si,  §i  are  the  singular  values  of  the  jth  px  p  block  of  the 
original  and  the  observed  image  matrices,  respectively.  In  [8], 
it  was  also  shown  that  the  Msvd  measure  gives  a  much  better 
representation  of  the  quality  of  an  image  with  regard  to  the 
HVS  for  a  large  variety  of  distortions.  Another  point  to  note 
is  that  a  lower  value  of  Msvd  indicates  better  image  visual 
quality.  For  the  same  experiment  which  generated  the  curves 
in  Fig.  2,  the  variation  of  the  Msvd  metric  for  output  image 
quality  is  as  shown  in  Fig.  3.  It  is  seen  that  as  the  threshold 


xl(r2 


Fig.  3.  Msvd  metric  based  assessment  of  SGWSR  HR  image  quality 


It  is  however  not  practical  (in  terms  of  computational  effort) 
to  plot  such  a  curve  every  time  the  SGW-based  (or  any 
other)  superresolution  algorithm  is  applied.  Furthermore,  in 
most  real-world  cases,  the  original  image  is  unavailable  for 
computing  the  Msvd  or  any  other  metric.  The  threshold  £ opt , 
which  achieves  the  optimal  trade-off  between  noise  reduction 
and  degradation  of  visual  quality  due  to  introduced  blur  is 
essentially  dependent  on  the  input  noise,  under  the  assumption 
of  similar  image  capture  devices  and  imaging  conditions2.  It 
follows  that,  for  a  specified  input  noise  level,  Copt.  determined 
for  a  small  number  of  training  images  will  be  near  optimal  for 
other  images  under  similar  imaging  conditions. 

Consider  a  real-valued  n\  x  n2  (without  loss  of  generality, 
assume  n\  >  n2)  image  matrix,  F  =  [/[?',  j]\  of  rank  R  ^  n2, 
with  singular  value  decomposition  (SVD) 


utfv 


£  0 
0  0  ’ 


E  =  diag(si, . . . ,  sr),  (12) 


where  si  ^  ^  Sr  >  0  are  the  singular  values  of  F. 

Consider  the  corruption  of  the  image  matrix  F  by  a  zero- 


mean  AWGN  matrix  N  to  get  F  = 


/M 


=  F  +  N.  Let 


§1  ^  ^  §r  >  0  be  the  singular  values  of  F,  which  is  of 

rank  R.  Then,  from  [13,  pp.  203,  Corollary  4.10],  it  follows 


that 


max{|§j  -  Sj|}  s£  ||iV||2,  i  =  l,2,...,R,  (13) 

i 


where  ||  •  ||2  is  the  matrix  2-norm.  A  generalization  of  the  above 
result  is  due  to  Mirsky  (holds  for  any  unitarily  invariant  norm) 
which  in  the  case  of  Frobenius  norm  (denoted  by  the  subscript 
Fr  for  the  norm)  becomes  [13,  pp.  205,  Corollary  4.13] 


i=l,2,...,R.  (14) 


is  increased  from  0,  there  is  initially  (until  approximately 
0.15)  a  significant  improvement  in  image  quality  due  to  noise 
reduction,  but  beyond  this,  as  the  effect  of  blur  introduced 
due  to  thresholding  starts  to  dominate,  the  quality  of  the 
reconstructed  image  suffers  with  increasing  threshold.  Thus, 
thresholding  involves  a  trade-off  between  noise  reduction  and 
degree  of  blur  introduced  in  the  reconstructed  HR  image,  and 
hence,  there  should  exist  a  threshold  value,  £opt,  the  location 
of  the  minimum  of  the  curve  in  Fig.  3,  for  which  this  trade-off 
is  optimal.  In  other  words,  the  Msvd  metric  is  sensitive  to 
visual  quality  degradation  due  to  blur. 

IV.  Optimal  Threshold 

The  optimal  threshold,  Copt-  for  the  SGWSR  algorithms 
presented  in  [1],  [2],  and  in  general,  for  any  denoising  require¬ 
ment,  is  the  one  which  minimizes  the  Msvd  metric  defined 
in  (11).  Since  (11)  involves  the  singular  values  of  a  matrix 
(which  cannot,  in  general,  be  expressed  in  terms  of  the  matrix 
elements),  it  is  not  possible  to  derive  a  closed-form  analytic 
expression  for  Copt-  The  optimal  threshold  is  thus  obtained  by 
plotting  the  curve  of  Msvd  values  of  the  reconstructed  HR 
image  for  various  values  of  £  and  then  finding  the  threshold, 
Copt,  corresponding  to  the  location  of  its  minimum. 


Consequently,  from  (11), 

Dj  <  \\N\\Fr.  (15) 

Since  (11)  involves  only  Dj’s,  the  Msvd  metric  is  essentially 
bounded  by  the  norm  of  the  noise  matrix,  N  =  [h[i,  j]}. 
Further,  if  the  singular  values  of  F  are  distinct  (simple),  then 
under  the  assumptions  of  sufficiently  small  noise  [13,  pp.  264- 
266], 

Si  =  Si  +  ujNvi  +  0(\\N\\l),  (16) 

where  U  =  [u\  . . .  uni\  and  V  =  [iq  . . .  vn2\  in  the  SVD  in 
(12).  Ignoring  the  higher  order  terms,  the  expected  value  of 
the  sum  of  squares  of  the  errors  in  the  singular  values  is 

^{E(^-^)2}  =  ^{E(^)2}  =  ^  (1?) 

where  fi2  is  the  variance  of  the  noise.  Thus,  the  Msvd  metric 
and  consequently,  the  optimal  threshold,  £opt,  is  dependent, 
on  the  average,  only  on  the  imaging  device  and  imaging  con¬ 
ditions.  Therefore,  for  a  given  class  of  devices  and  relatively 

“In  the  superresolution  algorithms  presented  in  [1],  [2],  the  input  blur  is 
due  to  the  image  capture  device  and  is  assumed  constant.  A  different  input 
blur  thus  constitutes  a  different  imaging  condition. 
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fixed  imaging  conditions,  the  desired  £opt  for  the  SGWSR 
algorithm  [1],  [2]  can  be  approximated  by  a  threshold  £opt 
computed  using  a  small  number  of  training  images. 

Simulation  results  which  support  the  analysis  presented 
above  are  given  in  Table  I.  70  different  HR  images  were 
reconstructed  from  their  respective  sets  of  6  LR  frames  which 
were  subjected  to  similar  input  degradations.  The  variance 
of  the  optimal  threshold  (which  was  computed  for  each  of 
these  HR  images  by  plotting  a  curve  and  finding  the  threshold 
corresponding  to  its  minimum  as  explained  earlier)  for  these 
images  is  noted  to  be  very  small.  As  seen  from  Table  I,  for 
another  set  of  70  images  with  different  imaging  conditions, 
the  optimal  threshold,  justifiably  different,  exhibits  very  low 
variance,  consistent  with  the  analytical  result  given  above. 


Fig.  4.  (a)  sample  LR  image;  (b)  original  image;  (c),  (d)  HR  images  based  on 
optimal  threshold  from  PSNR  and  Msvd,  respectively;  (e),  (f)  HR  images 
in  (c)  and  (d),  respectively,  after  post-deblurring  of  input  blur 


Training 

Set 

Imaging 

Conditions 

Variance 
of  C opt 

Copt 

(median) 

70 

AWGN  of  variance  0.001 
(5  x  5)  input  Gaussian  blur 
with  variance  1.0 

7.44  x  10“4 

0.02 

70 

AWGN  of  variance  0.01 
(5  X  5)  input  Gaussian  blur 
with  variance  0.5 

1.13  x  10“d 

0.125 

TABLE  I 

Variance  of  ( opt  for  different  imaging  conditions;  approximate 

OPTIMAL  THRESHOLD,  ( opt 


for  quantifying  visual  image  quality  prompted  the  use  of  the 
Msvd  metric  to  determine  the  threshold  for  optimal  trade-off 
between  the  blur  introduced  by  thresholding  and  input  noise 
reduction  in  the  reconstructed  HR  image.  Supporting  analysis 
and  simulation  results  are  given.  Subsequently,  the  use 
of  near-optimal  thresholds  (which  depend  on  the  imaging 
conditions  rather  than  the  images)  is  proposed  for  improving 
the  visual  quality  of  HR  reconstructions  in  the  SGWSR 
algorithms. 

Acknowledgement:  The  authors  thank  Dr.  Lina  Karam  and 
the  reviewers  for  their  prompt  and  constructive  evaluation. 


Finally,  a  sample  LR  frame  and  HR  reconstructions,  based 
on  the  SGWSR  algorithm,  with  the  optimal  thresholds  ob¬ 
tained  from  the  PSNR  curve  in  Fig.  2  and  the  Msvd  curve 
in  Fig.  3  are  shown  in  Fig.  4.  The  LR  images  were  subjected 
to  a  5  x  5  Gaussian  PSF  blur  and  corrupted  by  zero-mean 
AWGN  of  variance  0.01.  The  HR  images  obtained  from  the 
SGWSR  algorithm  based  on  the  optimal  threshold  derived 
from  the  PSNR  curve  and  the  Msvd  curve  are  shown  in 
Fig.  4(c)  and  (d)  respectively.  The  application  of  single  frame 
post-deblurring  (to  offset  the  effect  of  input  blur  as  well  as 
blur  introduced  by  thresholding)  based  on  the  accelerated 
Richardson-Lucy  algorithm  due  to  Biggs  and  Andrews  [14] 
(30  iterations)  to  each  HR  image  in  Fig.  4(c)  and  (d)  yields 
the  results  in  Fig.  4(e)  and  (f),  respectively.  The  PSNR- 
based  optimal  threshold  gives  better  noise  filtering  but  poor 
deblurring  [see  Fig.  4(c)  and  (e)]  while  the  Msvd- based 
optimal  threshold  produces  a  trade-off  between  introduced  blur 
and  input  noise  filtering  that  produces  output  images  with 
better  visual  quality  [see  Fig.  4(d)  and  (f)].  For  the  HR  image 
in  Fig.  4(d),  the  optimal  threshold  £opt,  as  per  the  Msvd  curve 
plotted  for  the  image  was  0.14.  From  Table  I,  the  near-optimal 
threshold,  £ opt  for  the  imaging  conditions  for  this  image  (row 
2  in  the  table)  is  0.125,  which  is  very  close  to  £opt.  The  image 
quality  due  to  the  use  of  £ opt  was  found  to  be  almost  identical 
to  the  ones  in  Fig.  4(d)  and  (f)  generated  using  £ opt ,  supporting 
the  use  of  near-optimal  thresholds. 

V.  Conclusion 

It  has  been  shown  that  the  choice  of  threshold  is  important 
in  the  SGWSR  algorithm.  The  inadequacy  of  PSNR  and  Q 
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Abstract 

The  main  contribution  of  this  paper  is  the  introduction  of  a  framework  for  estimation  of  multiple 
unknown  blurs  as  well  as  their  respective  supports.  Specifically,  the  Biggs- Andrews  (B-A)  mul¬ 
tichannel  iterative  blind  deconvolution  (IBD)  algorithm  is  modified  to  include  the  blur  support 
estimation  module  and  the  asymmetry  factor  for  the  Richardson-Lucy  (R-L)  update  based  IBD 
algorithm  is  calculated.  A  computational  complexity  assessment  of  the  implemented  modified 
IBD  is  made.  Simulations  conducted  on  real-world  and  synthetic  images  confirm  the  importance 
of  accurate  support  estimation  in  the  blind  superresolution  problem. 
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1  Introduction 


Blind  image  deconvolution  involves  image  restoration  from  degraded  observation(s)  with  either 
unknown  or  partially  known  information  on  the  type  and  extent  of  the  blur(s).  A  plethora  of 
single  channel  methods  have  been  documented  in  [1],  [2],  Among  methods  that  generalize  from 
single  frame  (channel)  blind  deconvolution  to  the  multiframe  (multichannel)  case  is  the  iterative 
one  proposed  by  Biggs  and  Andrews  [3].  Multichannel  blind  methods  that  do  not  have  single 
channel  equivalents  (however,  similarities  may  exist  as  in  [4])  have  also  been  proposed  [4],  [5].  The 
subspace  technique  [4]  is  non-iterative  and  like  other  approaches,  gives  erroneous  results  for  low 
signal-to-noise  ratio  (SNR)  cases  in  comparison  to  iterative  methods.  The  main  advantage  of  the 
multichannel  approach  arises  because  the  single  channel  methods  are  ill-posed  and  ill-conditioned. 

A  multichannel  blur  imaging  system  model  may  result  when  image  acquisition  is  through  multi¬ 
ple  cameras  as  shown  in  Figure  1.  Other  possibilities  leading  to  the  same  model  is  a  consequence  of 
image  capture  through  multiple  focuses  of  a  single  camera  or  images  acquired  from  a  single  camera 
through  a  changing  medium.  Phase  diversity  permits  restoration  of  target  and  unknown  point 
spread  function  by  two-channel  imaging,  typically  in-focus  and  out-of-focus.  There  is  a  natural 
multichannel  counterpart  [6].  A  ground-based  telescope  with  a  CCD  or  CMOS  sensor  can  be  used 
to  acquire  images  of  the  spot  in  the  solar  photosphere.  These  images  are  blurred  by  atmospheric 
turbulence  and  refractive  index  fluctuation  of  the  air  caused  by  temperature  variations. 

The  generic  model  in  Figure  1  includes  multichannel  acquisition  of  an  image  sequence,  a  module 
for  image  registration  based  on  the  projective  model  of  camera  motion  parameter  estimation  [7], 
[8],  a  module  that  simultaneously  implements  superresolution  (using  second  generation  wavelets) 
and  noise  filtering  by  either  soft  or  hard  thresholding  [9],  [10],  and  a  module  for  blind  blur  pararne- 
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ter  with  support  estimation  that  ultimately  produces  the  deblurred,  noise  filtered,  superresolved 


image.  Finite  support  linear  shift  invariant  blurs  are  reasonable  to  assume.  The  coprimeness 


Cameras  with  Multiple  Blurred,  Multiple  Blurred 

different  PSFs  Noisy  LR  Images  HR  Images 


D|Sf 


Figure  1:  Multichannel  Blind  Superresolution  Model 

condition  [11]  on  the  zeros  of  the  transfer  functions  representing  the  blurs,  needed  for  restoration 
from  multiple  deconvolution  operators,  if  satisfied,  is  expected  to  produce  a  high  quality  deblurred, 
noise  filtered,  superresolved  image.  This  coprimeness  condition  is  assumed  to  hold  in  the  blind 
identification  of  multichannel  finitely  supported  blurs  in  perfect  two-dimensional  signal  restoration 
[12],  A  multichannel  high  resolution  blind  image  restoration  scheme  was  also  considered  in  [13]. 
A  parameterized  blur  identification  (when  the  blurring  process  is  known  only  to  within  a  set  of 
parameters)  and  resolution  enhancement  scheme  which  includes  restoration  as  a  special  case,  has 
been  reported  in  [14],  The  iterative  blind  deconvolution  (IBD)  method  originally  proposed  in  [15] 
was  extended  in  [3]  to  multiple  frames  using  the  popular  Richardson-Lucy  (R-L)  algorithm. 

The  focus  of  this  paper  is  on  generalizing  the  result  in  [3]  to  include  not  only  multiple  blur 
identification  but  also  support  estimation  of  blurs  that  was  assumed  in  [3]  to  be  either  known  a 
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priori  or  determined  by  trial  and  error.  In  addition,  the  paper  also  presents  a  derivation  for  the 
asymmetry  factor  (discussed  in  greater  depth  in  Section  3)  required  in  the  IBD  algorithm  instead 
of  employing  a  trial  and  error  approach  as  proposed  in  [3],  [16].  It  is  also  pointed  out  that  the 
estimation  of  several  unknown  parameters  for  the  problem  of  reconstruction  of  a  high-resolution 
(HR)  image  from  multiple  under-sampled,  shifted,  degraded  frames  with  subpixel  displacement 
errors  (for  example,  the  model  advanced  by  Bose  and  Boo  [17])  has  been  considered  recently  in 
[18]- 

In  Section  2,  the  multichannel  point  spread  function  with  support  estimation  is  discussed.  The 
computational  complexity  of  the  procedure  presented  is  analyzed  and  simulation  results  on  real- 
world  as  well  as  synthetic  data  are  presented.  In  Section  3,  the  asymmetry  factor  in  the  modified 
Biggs-Andrews  (B-A)  multiframe  IBD  is  calculated  for  the  case  of  interest  here,  i.e.  procedure 
presented  using  R-L  IBD  algorithm  instead  of  least-squares  as  previously  done  in  [16].  In  Section 
4,  conclusions  are  drawn  and  the  new  results  obtained  in  this  paper  are  summarized. 

2  Multichannel  PSF  with  Support  Estimation 

The  resulting  HR  images  after  the  superresolution  step  in  Figure  1  still  need  to  be  deblurred  for  the 
best  possible  output  image  quality.  But,  in  most  practical  cases,  the  blur  of  the  acquisition  system, 
characterized  by  the  point  spread  function  (PSF)  is  unknown.  In  such  cases,  some  of  the  approaches 
that  can  be  adopted  are:  the  Biggs-Andrews  (B-A)  algorithm  [3],  the  You-Kaveh  regularization 
based  algorithm  [19]  and  the  Kundur-Hatzinakos  recursive  inverse  filtering  algorithm  [20] .  Though 
the  Kundur-Hatzinakos  algorithm  has  the  same  computational  complexity  as  the  B-A  algorithm, 
it  makes  several  restrictive  assumptions  (for  instance,  the  existence  of  the  inverse  of  the  original 
PSF)  while  the  complexity  of  the  You-Kaveh  approach  is  very  high.  A  more  detailed  comparison 
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of  these  and  other  techniques  is  presented  in  [21]. 


Among  the  three  mentioned  algorithms,  the  B-A  algorithm  generalizes  directly  to  the  multi¬ 
channel  case.  Since,  in  the  superresolution  problem,  extra  information  in  multiple  high  resolution 
frames  is  available,  it  is  desirable  to  use  a  method  which  can  exploit  it  for  improved  image  quality. 
Hence  the  B-A  algorithm  is  most  suited  for  the  multichannel  case  under  consideration.  Further, 
the  R-L  iteration,  which  constitutes  the  core  of  the  B-A  algorithm,  is  a  form  of  the  expectation 
maximization  (EM)  algorithm  of  Dempster  et  al.  [22]  (as  shown  by  Shepp  and  Vardi  [23]),  which  is 
an  efficient  and  widely  used  image  restoration  algorithm.  The  observed  image  is  obtained  by  stack¬ 
ing  in  3-D  the  degraded  2-D  observed  images.  The  procedure  then  becomes  similar  to  optimization 
algorithms  where  one  estimates  two  sets  of  parameters  in  the  iteration  by  freezing  one  and  updating 
the  other.  A  comment  relevant  to  this  type  of  technique  (alternating  minimization)  is  that  proof  of 
convergence  is  usually  not  feasible  without  appropriate  conditioning,  even  though  the  convergence 
of  the  EM  algorithm  and  hence  the  R-L  iteration  is  well  documented  [22],  [24],  However,  extensive 
simulations  have  confirmed  under  varying  conditions  (high  observation  noise,  varying  blur  support 
and  initial  over-estimate)  that  no  difficulty  in  rapid  convergence  was  witnessed.  Instability  is  not 
an  issue  when  using  R-L  iterations  since  it  can  be  shown  [3]  that  successive  iterates  are  bounded 
from  below  by  zero  and  from  above  by  the  total  energy  of  the  initial  estimate  (which  is  usually  the 
observation)  which  is  finite. 

2.1  The  Biggs- Andrews  Iterative  Blind  Deconvolution  (IBD)  algorithm 

The  B-A  IBD  algorithm  is  an  alternating  variable  optimization  approach  to  estimate  both  the 
image  and  the  PSF.  It  is  flexible  in  the  sense  that  different  optimization  techniques  can  be  applied 
for  the  image  and  the  PSF  updates  respectively.  Further,  the  algorithm  takes  into  account  different 
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rates  of  convergence  of  the  image  and  the  PSF  by  allowing  a  different  number  of  iterations  for  each 


within  the  main  alternating  optimization  loop.  This  approach  is  termed  as  ‘asymmetric  IBD’.  The 


essence  of  the  B-A  algorithm  is  captured  in  the  following  figure  (R-L  update  equations  are  used  for 


both  image  and  PSF):  The  basic  R-L  update  iterations  are  given  by  the  equations 


initial  estimates, 
f0  and  h0 


Figure  2:  Biggs- Andrews  (B-A)  IBD  algorithm 


h-k+l  —  ^  '  /fc  )  (1) 

Ln  Jk  \hk®  Jk  ) 

A+1  =  ' /“+1  *  Gs®a)  (2) 

where  g  is  the  original  observation,  •,  ★  and  <g>  represent  point-wise  multiplication,  correlation  and 
convolution,  respectively,  at  every  pixel  location  [i,j\  £  D  C  M2,  where  0  is  the  support  of  the 
image.  Similarly,  C  R2  is  the  support  of  the  PSF,  hk.  Acceleration  and  noise-dampening  as 
described  in  [3]  is  employed  in  the  implementation  though  it  is  not  shown  in  the  iterations  described 
in  Equations  (1)  and  (2)  above.  Further  details  about  the  algorithm  can  be  found  in  [3],  [16].  The 
initial  estimate  of  the  image,  /o,  is  set  to  be  equal  to  the  observation,  g.  The  initial  estimate  of 
the  PSF,  ho,  can  be  set  to  anything  with  the  exception  of  zeros  (one  of  the  properties  of  the  R-L 
algorithm  is  that  zero  values  remain  unchanged  over  iterations).  Biggs  and  Andrews  [3]  suggest  the 
use  of  the  autocorrelation  of  the  observed  image  as  the  initial  estimate  for  the  PSF.  In  the  results 
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presented  here,  a  matrix  with  all  values  equal  and  which  sums  to  unity  over  its  support  is  used  as 


the  initial  estimate  for  the  PSF.  Termination  of  the  iterations  occurs  when  either  the  norm  of  the 
difference  between  successive  estimates  falls  below  a  certain  preset  value  (ex:  ICY3)  or  a  maximum 
number  of  iterations  (ex:  30)  is  reached. 

Drawback:  The  biggest  drawback  of  the  above  mentioned  algorithm  is  that  it  requires  exact 
knowledge  of  the  support  of  the  PSF  for  optimal  performance.  Thus,  in  a  sense,  it  is  not  totally 
blind.  Typically,  since  there  is  no  information  about  the  PSF,  its  support  is  also  unknown.  In  [3], 
it  is  stated  that  when  this  information  is  unavailable,  the  only  alternative  is  to  over-estimate  the 
support  of  the  PSF  (so  that  it  contains  the  actual  support)  in  order  to  implement  the  deconvolution 
to  get  the  solution.  This  approach  does  not  yield  very  good  results  as  shown  by  simulation  results 
in  Section  2.3.  In  addition,  it  increases  the  computational  complexity  of  the  algorithm.  Under¬ 
estimation  of  the  PSF  support  leads  to  very  poor  deblurring. 

One  of  the  important  features  in  the  You-Kaveh  algorithm  is  that  it  includes  a  module  to 
iteratively  achieve  better  estimates  of  the  actual  support  of  the  PSF,  even  though,  like  the  B-A 
algorithm,  it  also  starts  with  an  over-estimation  of  the  support  of  the  PSF.  In  order  to  ensure  ease 
of  practical  implementation,  the  approach  in  [19]  involves  pruning  of  the  PSF  support  such  that  it 
is  always  rectangular.  Pruning  of  a  boundary-constituting  side  5(k  of  the  rectangular  PSF  support, 
Cfc,  is  carried  out  if 

hk{x)  ^  T,  V  x  €  5(k  (3) 

where  T  >  0  is  a  threshold  (the  threshold  should  be  positive  since  a  positivity  constraint  is  imposed 
on  both  the  image  and  the  PSF).  This  approach,  though  simple,  is  effective,  provided  a  proper 
threshold  is  chosen.  Various  simulations  presented  later  in  this  document  and  in  [19]  have  shown 
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that  this  method  reliably  estimates  the  support  of  the  PSF  after  a  small  number  of  iterations. 


2.2  Enhanced  Biggs-Andrews  IBD  with  PSF  Support  Estimation 

Since  the  You-Kaveh  algorithm  is  basically  designed  for  single  channel  deconvolution  and  better 
image  restoration  can  be  achieved  in  the  superresolution  problem  illustrated  in  Figure  1  by  making 
use  of  the  multiple  HR  images,  it  is  desirable  to  incorporate  the  iterative  PSF  support  estimation 
module  into  the  B-A  algorithm.  The  single  channel  case  is  first  considered  and  the  generalization 
to  the  multichannel  case  is  subsequently  presented. 

2.2.1  Single  Channel  PSF  Support  Estimation 

Since  the  support  estimation  needs  to  be  performed  only  on  the  PSF  (the  support  of  the  image 
is  known),  the  module  can  be  inserted  at  either  location  “1”  or  location  “2”  marked  in  Figure  2 
that  was  used  to  describe  the  B-A  algorithm.  The  following  arguments  present  the  advantages  of 
inserting  the  support  pruning  module  at  location  “2”  instead  of  location  “1” : 

•  Consistency  in  the  implementation  of  Equation  1  -  the  equation  assumes  that  fj~  is  the  image 
estimate  when  the  PSF  estimate  is  which  is  violated  to  a  certain  degree  if  is  pruned 
with  fk  remaining  unchanged.  On  the  other  hand,  insertion  of  the  module  at  location  “2” 
does  not  violate  any  of  the  equations.  See  Figure  3  for  the  modification  of  Figure  2  introduced 
here  so  that  the  blur  support  can  be  reliably  estimated.  Pruning  is  needed  because  of  the 
initial  overestimate  always  made. 

•  Initial  over-estimate  of  the  PSF  support  can,  possibly,  be  reduced  by  a  certain  degree  before 
an  image  update  is  performed 

The  updated  depiction  of  the  B-A  algorithm  is  as  shown  in  Figure  3  below.  The  condition  for 


initial  estimates, 
f0  and  ho 


Figure  3:  Modified  form  of  the  Biggs-Andrews  (B-A)  IBD  algorithm 

support  pruning  employed  in  the  algorithm  documented  in  this  article  is 

k 

^2  hj  (x)  ^  T,  X  <E  5(k  (4) 

j=k—t 

which  is  similar  to  that  employed  in  [19].  Here,  hk  represents  the  current  (kth  iteration)  PSF 
estimate,  T  the  threshold,  and  8£k  the  outermost  boundary  of  the  current  PSF  support,  Ck-  The 
main  differences  with  the  approach  in  [19]  are  that  the  amount  of  support  pruning  possible  in  each 
iteration  is  restricted,  and  an  option  to  include  information  about  variation  over  the  previous  t 
iterations  is  provided.  This  results  in  the  system  being  more  robust  to  temporary  fluctuations  in 
the  PSF  estimate.  The  threshold,  T,  is  heuristically  set  to  10%  of  the  energy  of  the  PSF.  Thus, 
though  the  rule  for  the  choice  of  the  threshold  is  determined  heuristically  and  is  independent  of  the 
type  and  support  of  the  PSF,  the  actual  value  of  the  threshold  depends  on  the  energy  of  the  PSF. 

2.2.2  Multichannel  PSF  Support  Estimation 

The  single  PSF  case  of  the  B-A  IBD  algorithm  can  be  easily  extended  to  the  multiple  PSF  case. 
That  is,  we  consider  imaging  process  given  by 

gn  =  f<S>hn,  n  =  1,2, . . . ,  N  (5) 


9 


where  /  is  the  original  image,  gn  represents  the  nth  blurred  observation  and  hn  is  the  corresponding 
PSF.  To  formulate  this  as  a  3-D  extension  of  the  single  channel  case,  it  should  fit  into  the  framework 
of  the  R-L  equations  ((1)  and  (2))  which  are  used  for  alternatively  updating  the  image  and  PSF 
estimates.  Essentially,  the  R-L  equation  is  an  iterative  maximum  likelihood  (ML)  approach  to 
estimate  one  of  the  terms  in  the  RHS  of 

g  =  f®h  (6) 

given  the  LHS  and  the  other  term  in  the  RHS.  Hence,  to  extend  the  IBD  algorithm  to  the  multi¬ 
channel  case,  Equation  (5)  needs  to  be  written  in  3-D  matrix  form  and  3-D  versions  of  the  operators 
in  Equations  (1)  and  (2)  are  to  be  used.  This  is  easily  achieved  as  follows: 

9[hj,n]  =  f[i,j,n\  <8>  h[i,j,n],  n  =  1,2, . . .  ,N  (7) 

where  <8>  represents  3-D  convolution.  The  values  of  g  and  h  for  fixed  values  of  the  3rd  dimension,  n, 
represent  the  N  blurred  observations  and  PSFs  respectively;  in  other  words,  the  blurred  observa¬ 
tions  and  PSFs  are  stacked  one  behind  the  other.  The  signal  /  can  now  be  viewed  as  a  3-D  signal 
with  f[i,j,  1]  being  the  original  image  and  with  f[i,j,  n]  =  0,  n  =  2, 3, . . . ,  N.  Thus,  the  R-L  equa¬ 
tion  can  be  applied,  and  consequently  IBD  can  be  carried  out  for  the  multichannel  case.  The  main 
difference  to  note  is  that  in  this  case,  the  initial  estimate  for  the  image,  /o,  and  the  observation, 
g,  are  not  the  same  since  the  initial  estimate  has  just  one  non-zero  frame  (to  maintain  the  validity 
of  Equation  (7)).  This  is  in  contrast  to  the  single  channel  case  where  the  initial  image  estimate  is, 
with  very  few  exceptions,  set  equal  to  the  observed  image.  The  initial  estimates  of  the  PSF’s  how¬ 
ever  are  set  to  be  matrices  with  equal  elements  summing  to  unity,  similar  to  the  single-channel  case. 

The  PSF  support  estimation  module  works  independently  on  each  of  the  PSF’s  and  iteratively 
arrives  at  the  correct  support  for  each  of  the  PSF’s,  though  different  number  of  iterations  may 
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be  required  for  convergence  to  the  true  support.  Further,  the  threshold  employed  for  support 
estimation  will,  in  general,  be  different  for  each  channel  since  the  threshold  is  set  to  be  10%  of  the 
energy  of  the  PSF.  The  next  subsection  presents  simulation  results  which  show  the  improvement  in 
restoration  quality  achieved  by  the  B-A  algorithm  with  the  use  of  the  support  estimation  module. 
It  also  demonstrates  the  accuracy  of  the  module  in  estimating  the  support  of  the  PSF. 

2.3  PSF  Support  Estimation:  Simulation  Results 

Results  for  the  single  channel  case  are  first  presented  followed  by  those  for  the  multichannel  case. 
The  results  show  the  importance  and  effectiveness  of  the  iterative  support  estimation  module.  Sin¬ 
gle  channel  results  are  presented  for  two  images:  one  real-world  image  and  the  other  a  synthetic 
image.  Simulation  results  for  the  synthetic  image  for  a  two  PSF  case  and  a  four  PSF  case  are 
subsequently  given.  The  synthetic  image  is  designed  such  that  it  presents  difficulties  to  the  restora¬ 
tion  algorithm  (due  to  the  very  sharp  contrast  and  transitions  in  the  image).  In  a  sense,  it  can  be 
thought  of  as  a  stress-test  for  the  algorithm.  Also,  different  PSF’s  are  used  for  generality. 

For  the  results  in  Figures  4(b)  and  4(d),  the  original  image,  PSF  (support  5x5)  and  the  blurred 
observation  are  shown  in  Figure  4(a).  The  results  in  Figure  4(d)  are  generated  using  the  original 
Biggs- Andrews  algorithm  with  different  guesses  for  the  support  of  the  PSF.  As  seen  from  the  figure, 
over  (9x9  and  7x7)  and  under  (3  x  3)  estimation  of  the  PSF  support  results  in  poor  quality 
of  restoration.  Over-estimation  results  in  inaccurate  restorations  -  the  grey  squares  in  the  original 
image  for  example,  appear  as  white  in  the  restoration.  Under-estimation  of  the  support  size  results 
in  the  output  remaining  significantly  blurred.  On  the  other  hand,  if  the  support  of  the  PSF  is 
guessed  correctly /known  a  priori  (5  x  5),  then  the  original  Biggs- Andrews  algorithm  yields  good 
restoration  results.  But  the  support  of  the  PSF  is  seldom  known  in  blind  deconvolution  problems 
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and  arriving  at  it  by  trial  and  error  is  impractical.  Figure  4(b)  shows  the  result  generated  for 
the  same  input  using  the  proposed  support  pruning  module.  Though  the  support  is  initially  over¬ 
estimated  as  9  x  9,  the  algorithm  iteratively  arrives  at  the  correct  support  as  shown  in  Figure  4(c). 
It  is  also  evident  that  the  restored  image  and  PSF  are  of  very  similar  quality  as  that  produced  by 
the  original  Biggs- Andrews  algorithm  when  the  support  is  known  a  priori  (5x5  case  in  Figure  4(d)). 

Connnents/inferences  similar  to  those  made  above  are  applicable  to  the  results  presented  for  the 
real-world  image  in  Figures  5(a),  5(b),  5(c)  and  5(d).  It  is  noted  that  over-estimation  of  the  support 
of  the  PSF  in  the  original  B-A  algorithm  (13  x  13  and  9x9  cases  in  Figure  5(d))  causes  excessive 
ringing  artifacts  at  the  boundaries  and  edges  while  under-estimation  (5x5  case  of  Figure  5(d)) 
results  in  insufficient  deblurring.  If  the  support  is  known  a  priori  or  arrived  at  by  trial  and  error, 
then  the  best  possible  quality  of  restoration  is  achieved  as  seen  from  the  7x7  case  in  Figure  5(d). 
On  the  other  hand,  if  the  proposed  algorithm  with  support  estimation  is  employed,  the  restoration 
is  of  good  quality  (Figure  5(b))  even  if  the  support  of  the  PSF  is  initially  over-estimated,  since  the 
algorithm  iteratively  arrives  at  the  correct  support  as  shown  in  Figure  5(c). 

In  the  multichannel  cases  presented  subsequently,  it  is  observed  that  the  PSF  support  estimation 
module  accurately  estimates  the  support  of  the  PSF  in  each  of  the  channels,  though  the  type  and 
support  of  the  PSF  in  each  channel  is  different.  Further,  it  is  also  noted  that  the  convergence  to 
the  true  support  in  each  channel,  as  expected,  takes  different  number  of  iterations.  The  quality  of 
the  restoration  is  also  much  better  than  in  the  single  channel  case  due  to  the  availability  of  more 
information  in  the  form  of  multiple  blurred  observations. 
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Synthetic  Image 
True  PSF 
Noise 


checkerboard  (64  x  64  pixels) 

5x5  Gaussian  with  variance  10 

additive  white  Gaussian,  zero  mean,  variance  0.0001 


Figure  4(a):  Original  image,  original  PSF,  and  noisy,  blurred  image  respectively 


Figure  4(b):  Restored  image  and  PSF  with  support  estimation  respectively 


Figure  4(c):  Iterative  PSF  support  estimation  (initial  over-estimate:  9x9) 
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Figure  4(d):  Restorations  with  different  guesses  for  PSF  support  (without  support  estimation) 
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Real-world  Image 
True  PSF 
Noise 


bridge  (512  x  512  pixels) 

7x7  Out-of-Focus/Circular 

additive  white  Gaussian,  zero  mean,  variance  0.0001 


Figure  5(a):  Original  image,  original  PSF,  and  noisy,  blurred  image  respectively 
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Figure  5(b):  Restored  image  and  PSF  with  support  estimation  respectively 


Figure  5(c):  Iterative  PSF  support  estimation  (initial  over-estimate:  13  x  13) 
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Figure  5(d):  Restorations  with  different  guesses  for  PSF  support  (without  support  estimation) 
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Synthetic  Image 
True  PSFs 

Noise 


checkerboard  (64  x  64  pixels) 

3x3  Gaussian  with  variance  10 
5x5  Out-of-Focus/Circular 

additive  white  Gaussian,  zero  mean,  variance  0.001 


Figure  6(b):  Restored  image  and  PSFs  with  support  estimation  respectively 
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Figure  6(c):  Iterative  support  estimates  for  the  two  restored  PSFs 
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Synthetic  Image 
True  PSFs 


Noise 


checkerboard  (64  x  64  pixels) 

7x7  Gaussian  with  variance  10 
5x5  Out-of-Focus/Circular 
7x7  Gaussian  with  variance  2.5 
3x3  Averaging 

additive  white  Gaussian,  zero  mean,  variance  0.001 


Figure  7(a):  Original  image,  original  PSFs  and  corresponding  noisy,  blurred  images  respectively 
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Figure  7(b):  Restored  image  and  PSFs  with  support  estimation  respectively 
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Figure  7(c):  Iterative  support  estimates  for  the  four  restored  PSFs 
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2.4  Complexity  Savings  due  to  Iterative  Support  Estimation 

It  has  been  shown  above  that  the  quality  of  the  deconvolved  output  in  the  B-A  algorithm  is  sensi¬ 
tive  to  the  support  of  the  initial  guess  for  the  PSF.  The  advantages  of  employing  an  iterative  PSF 
support  estimation  (by  pruning)  similar  to  the  You-Kaveh  algorithm  have  also  been  demonstrated. 
In  the  following  section,  the  savings  in  computation  due  to  the  iterative  PSF  support  estimation 
module  are  analyzed. 


The  basic  R-L  iterations  are  given  by  Equations  (1)  and  (2).  For  minimum  computational 
complexity,  the  convolution  and  correlation  operations  are  carried  out  in  the  Fourier  domain.  The 
above  equations  can  then  be  rewritten  as: 


hk+i  = 
fk+ 1  = 


1 


fk 

1 

hk+1 


hk  •  [Ft-T 


?-\Hk-Fk 


-fk  ■  T 


-l 


m+ 1-r 


F-\Hk+1  ■  Fk 


(8) 

(9) 


where  *  represents  complex  conjugation,  J-  and  J-  1  represent  the  Fourier  and  inverse  Fourier 
transform  operators  respectively,  and  X  denotes  the  Fourier  transform  of  the  signal  x. 


For  a  real  signal  of  length  N,  it  can  be  easily  shown  that  the  computation  of  the  Fourier  trans¬ 
form  and  its  inverse  require  0(N  log2(-/V))  operations.  It  can  also  be  shown  that  the  complexity  of 
an  n-dimensional  Fourier  transform  and  its  inverse  is  the  same  as  that  of  a  one-dimensional  signal 
of  length  equal  to  the  product  of  the  dimensions  of  the  n-dimensional  signal. 


Assuming  the  image  to  be  of  size  M  x  N  and  the  PSF  estimate  to  be  of  size  P  x  Q,  the 
complexity  of  implementing  each  of  the  terms  in  Equation  (8)  are  listed  below  based  on  the  above 
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paragraph  and  the  definitions 


C(X)^Cm(X)  +  Ca(X) 

Cm(X)  4  ^Xlog2(A)-5A  +  8 
Ca(X)  4  ^Xlog2(A)-5A  +  8 

where  Cm(X)  indicates  the  number  of  multiplications  required  for  an  A-point  FFT  and  Ca{X)  is 
the  corresponding  number  of  additions. 


List  of  Computation  Complexities 


computation  of  iq. 
computation  of  Hk 

(note  that  the  DFT  of  hk  should 
be  the  same  dimension  as  that  of  fk ) 
computation  of  tm\  =  T  1  (Hk  ■  Fk) 
computation  of  fm2  =  g/tm\ 
computation  of  tm^  =  T(tm‘}) 
computation  of  tm 4  =  J-~]  ( Ft *  ■  tm 3) 
computation  of  hk+i 


C(MN) 

C(MN) 


C(MN)  +  MN 

MN 

C(MN) 

C(MN)  +  MN 
2  PQ  +  MN 


The  list  for  implementing  Equation  (9)  is  the  same  as  above  except  for  the  last  item  on  the 
list,  whose  complexity  is  2 MN  +  PQ.  Thus,  the  PSF  and  image  update  equations  above  have  net 
complexities  0(5C(MN )  +  AMN  +  2 PQ)  and  0(5C(MN)  +  5 MN  +  PQ),  respectively.  Further, 
in  one  iteration  of  the  B-A  algorithm,  the  PSF  is  updated  n  times  and  the  image  m  times.  Let  K 
represent  the  total  number  of  iterations  of  the  B-A  algorithm.  Suppose  that  the  support  estimation 
routine  converges  to  the  true  size  of  the  PSF,  P  x  Q,  in  L  iterations  (for  simplicity,  it  is  assumed 


26 


that  there  is  a  dimensionality  reduction  of  2  in  each  dimension  for  each  of  the  L  iterations).  Then, 
the  savings  in  computational  complexity  due  to  the  support  estimation  routine  is 

O  (V  +  m)  ({K  -  L)(PQ  -  PQ)  +  +  Q))  ) 

For  the  multichannel  case,  the  savings  in  complexity  due  to  the  iterative  support  estimation  deter¬ 
mined  above  applies  to  all  the  PSF’s  (the  extent  of  savings  will  naturally  be  different  depending 
on  the  actual  support  of  the  individual  PSF’s). 

3  Asymmetry  Factor  for  R-L  IBD  Algorithm 

The  asymmetry  introduced  into  the  IBD  algorithm  by  Biggs  is  an  important  factor  deciding  the 
quality  of  the  deblurred  output  image  [16].  Biggs  (in  [16])  derived  an  expression  for  the  asymmetry 
factor,  AFk,  for  the  case  of  the  least-squares  error  metric  based  optimization  approach  to  blind 
deconvolution.  The  derived  expression  was  based  on  the  approximation  of  the  rate  of  convergence 
of  the  image  and  PSF  estimates  to  the  traces  of  the  corresponding  Hessian  matrices  for  the  error 
metric 

e  =  \\g  -  hk  ®  /fell2  (10) 

The  remainder  of  this  section  presents  a  method  of  computing  the  asymmetry  factor,  AF (.,  for  the 
case  of  the  R-L  IBD  as  an  alternative  to  employing  a  trial  and  error  based  approach  as  suggested 
in  [3],  [16]. 


The  metric  employed  here  is  the  log-likelihood  which  is  given  by 


£(g\fk,hk)  =  ^sTogrfc  -  rk 


(11) 


where  rk  is  the  reblurred  image  at  the  kth  iteration  used  in  the  R-L  update  equations  and  is  given 
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by 


rk  =  hk  <g>  fk  (12) 

For  notational  convenience,  the  Equation  (11)  is  written  as 

C  =  ^g{x)\ogrk(x)-rk{x)  (13) 

where  x  is  the  2-D  coordinate  [xq,  X2]  and  fl  C  M2  is  the  support  of  the  image,  i.e.  of  all  of  g,fk, 
and  rk.  Also, 

Tk{x)  =  hk(x )  <g>  fk{x) 

=  hk(t)fk{x  —  t)  (14) 

i£Cfc 

where  (k  C  M2  is  the  support  of  the  PSF  hk  and  t  is  the  2-D  coordinate  [ti ,  <2]  -  Further,  x  and  t  in 
the  above  expressions  are  assumed  to  have  a  lexicographical  ordering.  The  aim  now  is  to  construct 
the  Hessian  matrices  of  the  objective  function  given  by  Equation  (11)  with  respect  to  the  image 
and  the  PSF  estimate.  Consider  the  case  of  the  image  estimate  first.  The  partial  derivative  of  C 
with  respect  to  fk(y)  where  y  6  H  C  M2  is  given  by 

am  =  am{ 

(noting  that  g(x),  the  stack  of  observed  images,  is  independent  of  fk(y )) 
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Using  the  definition  of  rk(x)  from  Equation  (14) 


bMH)  <r‘(x,)  =  bMv)  1 5  ~  *> 


k 


d 


9fk{y) 


X  fk(t)hk(x-t) 


=  < 


(x-t)eCfc 
hk(x-y)  for  (pc  —  y)  6  Cfc 
0  otherwise 


(16) 


Hence,  substituting  Equation  (16)  in  Equation  (15), 


Apr) 

dfkiy)  ^  \rk(x) 
JkKyj  (x-y)G(k  v  v  7 


E 


-  1  hk(x  -  y) 


(17) 


In  the  above  derivation,  the  periodic  boundary  condition  was  assumed  for  the  image.  Proceeding 
further,  the  second  order  partial  derivative  with  respect  to  fk(z),  z  6  H,  is 


d2C 


d 


g{x) 

dfk(z)dfk(y)  dfk(z)  \  (t.J^ea  \rk(x) 


E 


-  1  hk(x  -  y) 


E  hk(X  -  y)-J— ( M  -  l 

(x-y)eCfc 

X]  hk(x-y)g{x) 

(x-y)€(k 


dfk(z)  \rk(x) 
d  (  1 


dfk(z)  \rk(x) 


=  ^-ht(x-Mx)W)sW)(n(x)) 

Using  the  result  for  the  partial  derivative  of  rk{ x)  obtained  in  Equation  (16), 


(18) 


d2c 


dfk{z)dfk(y) 


E 


-hk(x  -  y)hk(x  -  z)-^p- 

rk[x) 


(19) 


((x-y)£(k)n((x-z)£(k) 

The  above  equation  defines  the  Hessian  matrix  of  the  objective  function  given  by  Equation  (11)  with 
respect  to  the  image  estimate.  Due  to  the  presence  of  the  term  the  Hessian  matrix  generated 

by  the  generic  element  in  Equation  (19)  will,  in  general,  not  be  circulant  or  even  Toeplitz.  If 


=  1  (°r  a  constant),  then  the  matrix  is  a  circulant  [19].  Thus,  the  properties  of  circulant 
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matrices  cannot  be  exploited  (as  in  the  least-squares  error  case  considered  by  Biggs  in  [16])  for 
finding  the  trace  of  the  Hessian  matrix.  Only  the  diagonal  elements  need  to  be  computed  to 
find  the  trace  of  the  symmetric  matrix  generated  by  the  generic  element  in  Equation  (19).  This 
corresponds  to  the  case  z  =  y  in  equation  (19).  Therefore, 


d2C 


Y  - hk(x-y ) 


g(x) 

rKx ) 


(20) 


From  the  above  formula,  it  can  be  seen  that  the  RHS  of  Equation  (20)  will  contain  as  many  terms  as 
card(Cfc),  the  cardinality  of  the  support  of  the  PSF,  since  periodic  boundary  condition  is  assumed. 
Consequently,  the  trace  of  the  Hessian  matrix  will  be  given  by 


A/fc  -  Y 

yen 


d2C 

dW) 


Y  Y  ~hl(x  ~  y) 

yen  (: x-y)eCk 


g(x) 

r2(x) 


(21) 


The  expression  for  A fk  given  in  Equation  (21)  contains  all  the  elements  6  H.  Further,  with  a  little 
thought,  it  can  be  seen  from  Equation  (21)  that  each  element  of  H  will  appear  card(Oc)  times  with 
each  appearance  being  weighted  by  a  distinct  element  of  hk.  From  this  observation,  the  computa¬ 
tion  of  the  trace  is  relatively  straightforward. 


The  Hessian  matrix  of  the  objective  function  given  by  Equation  (11)  with  respect  to  the  PSF 
estimate  can  be  obtained  in  a  similar  manner  as  outlined  above  and  is  briefly  presented  below. 


DC 

dhk{y) 


d  ( 

777  ,  .  Y  9(x)  lo§  rk(x)  -  rk(x) 

aht (»)  vit 


E 


g{x) 


rk(x) 


-  1 


d 


dhk(y ) 


(xk(x)) 


(22) 


Using  the  definition  of  rk  (x)  from  Equation  (14) 


fk(x  -  y) 

0 


for  y  €  Cfc 
otherwise 


(23) 
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Hence,  substituting  Equation  (23)  in  Equation  (22) 


DC 

dhk{y ) 


1  fk(x  -  y) 


The  generic  element  of  the  Hessian  matrix  is  then  obtained  as 


(24) 


d2C 


dhk(z)dhk(y) 


Y  ~fk{x  -  y)g{ 


x)^r^  d ,  ,  (rk(x)) 


y&Ck 


rl{x )  d hk(z) 


Y  ~fk{x -y)fk(x-  z) 


(s/e4)n(ze4) 


g{x) 

rl(x ) 


(25) 


Again,  due  to  the  presence  of  the  term  9}?\ ,  the  Hessian  matrix  generated  by  the  generic  element 

rk\x) 

in  Equation  (25)  will,  in  general,  neither  be  circulant  nor  be  Toeplitz.  As  in  Equations  (20)  and 
(21),  the  trace,  Akk,  of  the  Hessian  can  be  computed. 


The  asymmetry  factor,  AFk  is  then  computed  as  given  by  Biggs: 


AFk  =  Eoo 


^4  hk  J_ 

fk  7 k 


where  7^ 


and  Eoo 


lA/fc  I 

^4  h2  Eo  / 

Eo  P  E<7  h 


(26) 


with  /  and  h  representing  the  original  image  and  PSF.  Since  this  information  is  unknown,  an 
estimate  of  the  value  of  E 00  is  usually  employed.  As  stated  by  Biggs  [16],  an  over  estimate  yields 
much  closer  performance  to  the  optimum  than  an  under  estimate.  The  R-L  IBD  algorithm  is 
pictorially  depicted  by  Figure  8  (refer  [16]). 


3.1  Asymmetry  Factor:  Simulation  Results 

A  simulation  result  is  presented  to  compare  the  performance  of  the  R-L  IBD  with  automatic  in-loop 
estimation  of  the  asymmetry  factor  AFk,  against  that  of  the  R-L  IBD  with  manual  tuning  of  the 
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initial  estimates  f0,  h0 
k=  0 


Figure  8:  Modified  R-L  update  based  Asymmetric  IBD  algorithm 

fixed  asymmetry  ratio.  The  four  channel  case  presented  earlier  is  considered  for  comparison.  Figure 
9  shows  the  restoration  results  for  both  cases,  from  which  it  is  clear  that  the  performance  of  the 
R-L  IBD  with  automatic  in-loop  estimation  of  the  asymmetry  factor  AF^  closely  matches  that  of 
the  optimally  tuned  R-L  IBD  with  a  fixed  asymmetry  ratio. 


4  Conclusions 

In  the  IBD  algorithm  shown  in  Figure  2,  the  integers  m  and  n  which  are  representative  of  the 
convergence  rates  of  the  image  and  point  spread  function  updates,  respectively,  are  not  equal 
(asymmetric  IBD)  and  the  image-dependent  asymmetry  factor  n/m  was  selected  by  trial-and-error 
(manual  tuning).  The  main  contributions  of  this  paper  are:  first,  the  insertion  of  a  support  estima¬ 
tion  (by  pruning)  module  for  each  finitely  supported  PSF  as  shown  in  Figure  3,  and,  second,  the 
computation  of  the  asymmetry  factor  by  finding  the  traces  of  the  respective  Hessian  matrices  whose 
generic  elements  were  calculated  in  Equations  (19)  and  (25).  The  computation  of  the  traces,  natu- 
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manually  tuned 


original 


automatic 
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rally,  involves  the  finding  of  only  the  diagonal  elements  of  these  matrices,  as  given  in  Equation  (20) 
and  its  counterpart  that  is  readily  derivable  from  Equation  (25).  These  contributions  to  multiframe 
blind  deconvolution  of  image  sequences  degraded  by  multiple  blurs  in  conjunction  with  the  already 
developed  algorithms  for  superresolution  and  simultaneous  noise  filtering  using  second  generation 
wavelets  and  thresholding  [9],  [10],  as  shown  in  Figure  1,  will  contribute  significantly  towards  the 
goal  of  obtaining  a  superresolved,  deblurred  and  noise  filtered  image  from  image  sequences  acquired 
with  video  cameras. 
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Superresolution  and  Noise  Filtering  Using  Moving 

Least  Squares 

N.  K.  Bose,  Life  Fellow,  IEEE  and  Nilesh  A.  Ahuja 


Abstract — An  irregularly  spaced  sampling  raster  formed  from 
a  sequence  of  low  resolution  frames  is  the  input  to  an  image 
sequence  superresolution  algorithm  whose  output  is  the  set  of 
image  intensity  values  at  the  desired  high  resolution  image  grid. 
The  method  of  moving  least  squares  (MLS)  in  polynomial  space 
has  proved  to  be  useful  in  filtering  the  noise  and  approximating 
scattered  data  by  minimizing  a  weighted  mean  square  error 
norm,  but  introducing  blur  in  the  process.  Starting  with  the 
continuous  version  of  the  MLS,  an  explicit  expression  for  the 
filter  bandwidth  is  obtained  as  a  function  of  the  polynomial 
order  of  approximation  and  the  standard  deviation  (scale)  of 
the  Gaussian  weight  function.  A  discrete  implementation  of  the 
MLS  is  performed  on  images  and  the  effect  of  choice  of  the  two 
dependent  parameters,  scale  and  order,  on  noise  filtering  and 
reduction  of  blur  introduced  during  the  MLS  process  is  studied. 

Index  Terms — moving  least  squares,  Hermite  polynomials, 
superresolution 

I.  Introduction 

Most  approaches  to  noise  filtering  in  images,  ranging  from 
weighted  least  squares  [1]  and  weighted  total  least  squares  [2] 
to  bilateral  filtering  [3]  assess  image  quality  by  calculating 
either  the  mean  square  error  (MSE)  or  the  peak  signal-to- 
noise  ratio  (PSNR).  Other  measures  of  visual  quality  have 
been  considered  including  most  recently  [4]  a  measure  based 
on  singular  value  decomposition  (Mg^rO  of  the  image  matrix, 
introduced  first  in  [5],  Noise  filtering  with  superresolution 
has  been  implemented  in  several  approaches  starting  with 
[1]  and  leading  upto  current  approaches  [6],  The  problem  on 
deblurring  the  input  blur  present  in  low-resolution  (LR)  frames 
has  been  tackled  in  various  ways  [2]  [7],  Post-processing 
has  been  the  main  recourse  for  input  blur  removal  in  blind 
superresolution,  where  the  blurs  have  to  be  estimated.  Note 
that  the  blur  parameters  and  the  respective  blur  supports  can 
be  estimated  from  the  LR  frames  in  a  pre-processing  step  [8, 
and  references  therein],  but  the  actual  deblurring  should  be 
performed  after  the  image  fusion  as  a  post-processing  step 
[7],  [9].  In  a  recent  paper,  the  idea  of  implementing  deblurring 
and  interpolation  simultaneously  has  been  pursued  [10]  in  the 
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case  of  known  or  a  priori  estimated  blurs  (atmospheric  and 
camera)  and  an  underdetermined  system  of  linear  equations 
that  result  from  minimization.  This  underdetermined  system 
occurs  when  the  number  of  frames  Nf  <  r2,  where  r  is 
the  resolution  enhancement  factor  along  each  direction.  Blurs 
can  also  appear  during  denoising  of  an  unblurred  but  noisy 
input.  For  example,  in  [4]  and  [6],  the  promising  technique 
of  wavelet-coefficient  thresholding  introduces  undesirable  blur 
necessitating  the  choice  of  an  optimal  threshold  (the  only 
parameter)  for  tradeoff  between  blur  removal  (introduced  by 
thresholding)  and  noise  filtering.  The  objective  of  this  paper 
is  to  continue  the  investigation  into  tradeoff  between  the 
opposing  effects  of  noise  and  introduced  blur  in  other  methods 
that  apply  to  multiframe  superresolution.  The  method  of  choice 
for  such  a  study  in  this  paper  is  that  of  (weighted)  moving  least 
squares  recently  used  in  single  frame  image  processing  [11] 
[12]. 

Moving  Least  Squares  (MLS)  is  a  local  approximation 
method  that  has  been  documented  in  [13,  Ch.4].  Here,  its  scope 
in  a  problem  of  current  interest,  namely  superresolution,  is 
investigated.  Consider  the  irregularly  sampled  raster  in  Figure 
1  generated  by  the  sequence  of  LR  frames  captured  by  a  video 
camera  whose  motion  parameters  may  be  described  by  the 
general  projective  model  [14].  The  goal  of  superresolution 
is  to  convert  this  irregularly  sampled  raster  in  Figure  1  to 
a  regularly  sampled  one  where  this  regularly  sampled  raster 
is  comprised  of  grid  points  at  equispaced  intervals  along 
each  of  the  coordinate  axis.  The  pixel  value  at  each  grid- 
point  is  computed  with  a  polynomial  approximant  using  the 
pixels  in  a  defined  neighborhood  of  the  grid-point  under 
consideration.  Since  the  defined  neighborhood  may  change 
from  one  grid-point  to  the  next,  the  coefficients  as  well  as  the 
order  of  the  polynomial  approximant  has  to  be  adaptive.  The 
MLS  approach  to  superresolution,  discussed  here,  mandates 
the  choice  of  two  parameters  for  optimal  trade-off  between 
introduced  blur  and  noise  filtering.  The  two  parameters  are  the 
order  of  approximation  and  a  scale  parameter  characterizing 
the  standard  deviation  of  the  Gaussian  weight  function  for 
the  predictor  function  coefficients  as  described  next.  For  an 
overview  of  various  superresolution  techniques,  see  [15], 


II.  Continuous  Moving  Least  Squares  As  An 
Equivalent  Filtering  Operation 

For  (/-variate  polynomial  interpolation,  the  space  of  interest 
is  the  space  II  ^  of  ri-variate  polynomials  of  total  degree  <  s. 
Since  the  dimension  of  this  polynomial  space  is  (for  a  proof. 
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•  ■  *  Samples  of  low-resolution  images 

Fig.  1.  Irregularly  sampled  raster  superimposed  on  HR  grid 

see  [13,  pp.20-21]) 


multivariate  interpolation  requires  the  consideration  of  a  set 
of  sample  points  (nodes)  of  cardinality  at  least  Mf.  The 
problem  of  characterizing  nodes  that  admit  unique  Lagrange 
interpolation  has  been  extensively  documented.  Hermite  in¬ 
terpolation  allows  coalescence  of  nodes  along  a  line  while 
coalescence  of  nodes  along  a  plane  gives  rise  to  the  more 
general  Birkhoff  interpolation.  It  is  easy  to  show  that,  unlike 
in  the  univariate  case,  a  set  of  distinct  points  may  not  admit 
a  unique  Lagrange  interpolation  in  the  multivariate  case. 
For  interpolation  by  multivariate  polynomials  in  Ilf  to  be 
possible,  the  data  sets  have  to  be  restricted  [13,  pp.  19] .  Due 
to  the  fundamental  complexities  in  multivariate  interpolation, 
it  is  expedient  to  sacrifice  the  global  approach  in  favor  of 
local  surface  approximation  approaches  as  was  done  in  [7] 
and  the  MLS  approach  [11]  [12]  for  image  approximation 
and  interpolation  (bivariate  case).  In  the  MLS  method,  the 
optimal  fitting  function  is  expressed  as  a  linear  combination 
of  basis  functions,  whose  coefficients  are  chosen  to  minimize 
the  weighted  mean  squared  error  (MSE)  between  the  signal 
and  its  approximant.  The  MSE  is  generated  by  moving  the 
fitting  function  coefficient  mask  from  position  to  position  in 
the  image  plane.  The  primary  difference  between  the  approach 
in  [11]  [12]  and  the  approach  in  this  paper  is  the  incorporation 
of  a  variable  scale  in  the  weight  function  and  variable  order 
of  approximant  so  as  to  capture  better  the  distribution  of 
irregularly  spaced  samples  in  the  problem. 

To  analyze  the  effect  of  the  choice  of  scale  and  order  on 
blur  and  noise  in  the  reconstructed  high-resolution  image,  a 
continuous  formulation  of  MLS  is  used.  Moreover,  for  the 
sake  of  simplicity,  the  analysis  is  initially  carried  out  for  a 
1-D  signal.  The  generalization  to  2-D,  in  product  separable 
form,  is  reasonably  straightforward  and  is  explained  later  in 
this  section. 

A.  Analysis  of  1  -D  Case 

Here,  we  summarize,  first,  the  results  described  by  Boom- 
gard  and  Weijer  [11]  and  Fenn  and  Steidl  [12],  before  deducing 
the  formula  for  filter  bandwidth  as  a  function  of  the  two 
parameters,  scale  and  order.  In  the  d  =  1  case,  let  the  function 


to  be  approximated  by  MLS  be  f(x).  Locally,  about  a  fixed 
but  arbitrary  point  xo,  the  function  is  approximated  as  a  linear 
combination  of  basis  functions  (f>o(x).  cf>i(x), . . . ,  <j>k{x).  If  the 
basis  functions  are  the  monomials  l,x,  ...,xk  respectively, 
then  this  approximation  can  be  viewed  as  the  truncated  Taylor 
series  expansion  of  f(x)  about  x  xq.  The  approximation  of 
f(x)  about  x  =  Xq  is  denoted  as  f(x).  Then,  the  fcth  order 
approximant  is, 

/( x)  =  a0(x0)(j)0(x  -  a;0)  H - b  ak{xo)(j>k{x  -  cc0).  (1) 


The  coefficients  in  the  above  expansion  depend  on  Xq.  For 
the  sake  of  brevity,  the  argument  xo  will  be  dropped  from 
the  notation  and  the  coefficients  will  be  referred  to  simply  as 
a,i,i  =  0 , ,k.  These  coefficients  are  found  by  minimizing 
the  weighted  norm  \\f(x )  —  /(a;)||^,  which  is  the  weighted 
MSE.  The  weighting  function,  w(x),  is  such  that  the  points 
close  to  the  original  ,c0  are  weighted  higher  than  the  points 
far  from  x$.  Thus  the  approximation  of  /( x)  by  f(x)  is  made 
local  in  nature  by  the  choice  of  w(x).  A  typical  choice  for 
w(x)  (ubiquitous  in  adaptive  system  theory,  e.g.  [16])  is  the 
Gaussian  function 

w(x)  =  ^—e-*2/2°2,  (2) 

ctv27t 

having  scale  (standard  deviation)  a.  From  standard  least- 
squares  theory,  the  set  of  coefficients  a,i,i  =  0, . . . ,  k  that 
minimizes  the  weighted  MSE  is  the  solution  to  the  following 
system  of  linear  equations  [11]: 


(/;  — 
{4>i, 


I  f{u)(j)i{u  -  x0)w(u  -  X0)du 

J  —  OO 

/oo 

<j>i(u  —  Xo )4>j{u  —  Xo )w(u  —  Xo )du  . 

-OO 


Instead  of  the  monomials  1,  x, . . . ,  xk,  scaled  Hermite  poly¬ 
nomials  of  order  upto  k  can  be  chosen  as  the  basis  functions. 
The  Hermite  polynomial,  Hn  (®),  of  order  n  is  defined  as 

Hn{x)  =  {-  l)"e*2A(e-*2).  (4) 

Hermite  polynomials  upto  order  k  span  the  same  space  as 
monomials  upto  order  k,  but  are  orthogonal  with  respect  to 
the  Gaussian  weighting  function  w(x),  i.e. 


Hm(x)Hn(x)e  x  dx  = 


where  5mn  is  the  Kronecker  delta  function.  By  choosing  the 
basis  functions  as  4>k(x)  =  Hk(x/ay/2),  the  matrix  (called 
Gram  matrix)  in  Eq.  (3)  becomes  diagonal  [12]  and  the 
diagonal  elements  are  ||</>i||^  =  2lil.  Using  this  in  Eq.  (3), 
the  coefficients  oq,  *  =  0, . . . ,  k  are  easily  shown  to  be 


_  _  (/ ! 

Hi\\l  2?;*! 


(5) 
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Then,  from  Eq.  (1), 


f(x)  =  E  ai<^(x  ~  ^o)  =  e 


(/)  <f>i)i 
2*?:! 


fa(x  -  X0). 


i—0  i— 0 

The  equation  above  is  an  approximation  of  /(&),  locally 
around  x  =  Xq.  The  value  of  the  approximant  f(x o)  is 
obtained  by  evaluating  the  preceding  equation  at  x  =  Xq  to 
get 


f{x  o)  =  E 

Define 


i= o 


(f,  4>i)t 

2  H\ 


i<(o  )  =  (/,E 


2  =  0 


»i(o  )4>i 

2  H\ 


Kx )  =  E 


i—0 


fa(Q)fa(x) 

2H\ 


(6) 


Then, 


/(^o)  =  (/,  /i)w  = 


'  — OO 
pOO 


f(u)h(u  —  Xo  )w(u  —  Xq  )du 
f(x  o  +  u)h(u)w(u)du , 


ff(x)  =  h(x)w(x)  =  E 


i= 0 


<fc(0)(fc(x) 

2*i! 


w(x). 


A2  = 


1 

2a2 


1  + 


4EUEL 


EZ 

i—0  2-/  7=0 


i=0  Z-^j—0  8i+Ji!j!(i+j)! 


[2(»+j)]i 


(2):  i 

w(x,y)  =  - — Te  ^2+2/2)/2ct2  =  w(x)w(y).  (10) 

Z7raz 

The  basis  functions  faj(x,y)  are  chosen  as  the  bivariate 
Hermite  polynomials  HlJ  0/ crv/(2)>y/<J\/(2))  in  product 
separable  form, 

<j>ij(x,y)  =  fa(x)fa(y).  (11) 

The  sum  of  the  highest  degrees  of  fa  (x)  and  fa  ( y )  chosen  is 
the  order  of  approximation.  Here,  the  order  of  approximation 
is  2k.  The  value  of  the  function  at  x0  is  then  estimated  by 
evaluating  the  polynomial  function  in  Eq.  (9)  at  x  =  Xo. 
Similar  to  Eq.  (5),  the  coefficients  al3  can  be  found  as 

_  (/)  * Pij)w 


ll'/hj  II  w 

Using  Eq.  (12)  in  Eq.  (9)  results  in 

k  k 

f{x,y)  =  EE^-^') 

*= o  3=0 


(12) 


faj(x  -x0,y-  yo) 


2 

13  1 1  w 


and  f{x o)  can  be  obtained  by  linear  shift-invariant  filtering 
of  /( Xq  +  x)  with  a  system  characterized  by  the  unit  impulse 
response  h(x)w(x).  It  is  desired  to  find  the  bandwidth  of  the 
filter  whose  unit  impulse  response  is  (using  Eq.  6) 


which,  evaluated  at  (x,y)  =  (xo,yo)  yields, 

f(xo,Vo)  =  EE</’ 

i—0  j—0  WvUw 


k  k 


=  /-EE 


(7) 


fa  j  (0,0)  fa 

i=0  j=0 


(13) 


In  the  Appendix,  it  is  proved  that  the  bandwidth  A2  is 
explicitly  computable  in  terms  of  the  standard  deviation  a  and 
the  order  k  =  21  or  k  =  21  +  1  as  given  next 


Again,  as  in  the  1-D  case  of  Eq.  (7),  the  equivalent  filter  unit 
impulse  response  is 


k  k 


i—0  j—0 


faj(010)faj(x,y) 


>(x,y).  (14) 


(8) 


=0  Z^jj— 0  &+ii\j\(i-\-j)\ 

Before  proceeding  to  understand  the  implications  of  the  above 
expression,  the  linear  shift-invariant  system  theory  based 
analysis  performed  in  this  subsection  is  extended  to  2-D 
signals. 


Using  the  separability  of  the  basis  functions  and  the  weight 
function,  the  above  expression  simplifies  to 

g(x  y)=  (y  MMM wrx)\  ( y  ) 

-  Pi  (a;) 52  (y),  (15) 


B.  Generalization  to  Separable  2-D 

Suppose  that  it  is  desired  to  approximate  the  function  /(x)  = 
f(x,y)  at  a  point  xo  =  (xo,j/o)-  The  function  /(x)  can  be 
approximated  locally  around  Xo  by  a  linear  combination  of 
the  basis  functions  faj(x)  =  faj(x,y),  i  =  and 

j  =  1, . . . ,  fe.The  estimate  f(x,  y)  about  (xo,  yo)  is,  therefore, 
expressed  as, 

k  k 

f(x > */)  =  EE aij<l>ij(x  ~xo,y~  yo),  (9) 

i—0  j—0 

where,  without  loss  of  generality,  ki  =  k2  =  k.  Like  in  the  1-D 
case,  the  coefficients  are  found  by  minimizing  the  weighted 
norm  \\f(x,y)  —  f(x,y)  ||^,  where  the  weighting  function, 
w(x,y),  is  chosen  to  be  product  separable  counterpart  of  Eq. 


where  g\(x)  and  g2(y)  are  the  summations  in  the  respective 
brackets.  The  2-D  Fourier  transform  (FT)  of  g(x,  y)  is 


/oo  /»oo 

/  gi(x)g2(y)e~:,(UlX+U2y)dxdy 

-OO  J  —  OO 


gYY^dx 


—  Gl(u>l)G2(u>2 ), 


g2(y)e  JU>2Vdy 

(16) 


where  G\(uj\)  and  G2(u>2)  are  the  integrals  in  the  respective 
brackets.  Finally,  the  radial  bandwidth  of  this  filter  is 

,2  _  fZofZoiv?  +  (xZ)\G(u1,u2)\2duj1duj2 

f-oof-oo  \G(u)1,U2)\2dLU1dLU2 

=  A2  +  A2 

iOl  '  UJ2  ’ 
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where, 

2  _  f^oof^oo^ l|Gi(cJi)|2|G2(^2)|2^1^2 

JZoIZo  |Gi(wi)|2|G2(o;2)|2<iwidw2 
_  f-o 0ujf\G1(uJi)\2duj1 
/_0000|G1(a;1)PdWl  ' 

and  A2,  ,  is  defined  likewise.  Thus,  the  1-D  result  is  sufficient 
to  evaluate  the  separable  2-D  case  too. 

C.  Inferences 

The  analysis  of  the  previous  two  subsections  shows  that 
the  MLS  technique  can  be  viewed  as  a  filtering  operation.  An 
expression  for  the  bandwidth  of  this  filter  is  given  in  Eq.  (8). 
The  following  observations  can  be  made  from  this  expression: 

1)  An  increase  in  scale  causes  a  decrease  in  bandwidth 
and  vice  versa.  This  result  is  intuitively  appealing  since 
a  larger  scale  means  that  the  support  of  the  weighting 
function  w(x)  and  hence  of  the  filter  unit  impulse 
response  g(x)  is  larger,  which  causes  greater  averaging 
or  smoothing  (blurring).  In  the  frequency  domain,  this  is 
equivalent  to  a  reduction  in  bandwidth,  which  produces 
more  noise  filtering. 

2)  The  expression  for  dependence  on  order  is  seemingly 
complicated.  A  plot  of  filter  bandwidth  at  a  fixed  scale 
versus  the  order,  described  in  Eq.  (8),  is  shown  in  Figure 
2.  The  plot  shows  that  bandwidth  increases  approxi¬ 
mately  linearly  with  filter  order.  This  is  not  surprising 
either  since  by  using  a  higher  order,  rapid  changes  (or 
high-frequency  components)  in  the  original  signal  can  be 
modelled  more  accurately.  Again,  in  frequency  domain 
terms,  this  is  equivalent  to  an  increase  in  bandwidth  that 
admits  more  noise. 

For  image  sequence  superresolution,  the  above  two  effects 
result  in  a  tradeoff  between  noise  in  the  reconstructed  image 
and  the  blur  introduced  by  the  MLS  process.  Increasing  scale 
lowers  the  bandwidth  of  the  filter  resulting  in  greater  noise 
filtering  but  more  output  blur.  Increasing  order  increases  the 
bandwidth  of  the  filter  thereby  reducing  the  blur  introduced 
by  the  process  at  the  expense  of  increased  noise  in  the 
reconstructed  image. 

III.  MLS  Applied  To  Image  Superresolution 

In  image  superresolution,  the  input  is  known  at  a  set 
of  irregularly  spaced  discrete  data  points  generated  from  a 
sequence  of  registered  LR  frames  and  the  aim  is  to  estimate 
the  image  intensity  on  the  regularly  sampled  raster  of  the  HR 
grid.  For  this,  a  discrete  implementation  of  the  MLS  procedure 
is  used.  In  this  section,  the  discrete  formulation  of  MLS 
is  explained,  its  similarity  to  local  polynomial  regression  in 
statistics  is  pointed  out  and  a  technique  to  adaptively  select  the 
scale  and  order  parameters  for  estimating  the  image  intensity 
at  each  point  of  the  HR  grid  is  described. 

A.  Discrete  MLS 

Given  the  samples  of  a  function  /(x)  =  f(x  1,2:2)  at 
N  points  xm  =  (xmi ,  xrn2 ) ,  m  =  1,...,  A,  it  is  de¬ 


Filter  bandwidth  vs  order/2, a  =1 


Fig.  2.  Filter  bandwidth  vs  order/2 


sired  to  estimate  the  value  of  the  function  at  some  ar¬ 
bitrary  but  fixed  point  X(j.  Locally  about  Xo,  the  func¬ 
tion  is  expressed  a  linear  combination  of  basis  functions 
(f> o(x),  (j>i  (x), . . . ,  </>fc(x).  If  the  basis  functions  are  the  bivari¬ 
ate  monomials  1,  X\,  X2,  xf,  X1X2,  x\, . . .  ,etc.,  then  this  can  be 
viewed  as  the  Taylor  series  expansion  of  /(x)  about  x  =  x(). 
Denote  the  approximation  of  /(x)  about  x  =  Xo  by  /(x). 
Then, 

/(x)  =  a04> 0(x  -  x0)  -I - 1-  afe^fc(x  -  x0).  (17) 

The  coefficients  a,i,  i  =  0, . . . ,  k  are  found  by  minimizing  the 
weighted  norm  ||/(x)  —  /(x)||^,  where  the  weight  function 
w(x)  =  w(xi,X2 )  is  chosen  to  be  the  circularly  symmetric 
bivariate  Gaussian  function  given  by 

W(X)  =  2^e~M2/2<T2’ 11x1,2  =  X*  +  Xl  (18) 
This  results  in  a  system  of  linear  equations  similar  to  the  one 
in  (3).  However,  the  weighted  inner-product  and  the  weighted 
norm  are  now  defined  differently  as 

N 

(< k,4>i)w  =  ^2  fi(Xm  -x0)</>j(xm  -x0)w(xm  -x0) 

m=  1 

=  (19) 

N 

\\<t>i\\w  =  I <MXm  -  x0)|2w(xm  -  X0) 

771=1 

=  (20) 

where, 

=  [<Mxi  -  xo)  •  ■  ■  </>i(xA  -  xo)]T  (21a) 

W  =  diag(w(xi  -  x0), . .  .,w(x N  -  x0)).  (21b) 

Thus,  the  counterpart  of  Eq.  (3)  for  discrete  MLS  can  be 
compactly  represented  as 

$TWf  =  $TW<ka, 


(22) 
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where. 


(  /(X  l)  \ 

(  a°  \ 

— c 

o 

II 

f  = 

^  /(xjv)  j 

,  a  = 

\  ak  J 

and  the  solution  of  the  system  is,  therefore,  given  by 

a=($TWf)'1$TWf,  (23) 

where  the  inverse  of  <l>7  W<1>  exists  under  mild  constraints  on 
the  number  N  and  distribution  of  the  irregularly  spaced  points 
about  xo-  The  above  technique  is  closely  related  to  that  of 
local  polynomial  regression  found  in  computational  statistics 
[17],  Local  regression  is  used  to  model  a  relation  between  a 
predictor  variable  x  and  a  response  variable  y  related  to  the 
predictor  variable.  Suppose  a  dataset  consists  of  A’-pairs  of 
observations,  (xi,  j/i), . . . ,  (x.n>  Un )•  A  model  of  the  form 

ym  =  Mxm)  +  fm  (24) 

is  assumed  where  //(x)  is  an  unknown  function  and  em  is  an 
error  term  representing  the  random  errors  in  the  observations. 
The  errors  em  are  assumed  to  be  independent  and  identically 
distributed,!. e.  the  covariance  matrix  of  the  errors  is  a2 I  where 
cr2  is  the  variance  of  each  of  the  error  terms.  Locally,  around 
a  point  xo,  it  is  assumed  that  p(x)  can  be  well  approximated 
by  /t(x)  in  terms  of  its  Taylor  series  expansion  up  to  a  certain 
order  k,  i.e. 

A(x)  =  a0<Mx  -  xo)  H - h  ak(j)k(x  -  x0), 

where  </>o(x),  ^i(x),  02  (x),  ■  ■  •  are  the  bivariate  monomials 
1,  xi,  x-2, . . .  respectively.  By  following  a  weighted  least 
squares  technique  similar  to  the  one  used  for  discrete  MLS, 
an  expression  similar  to  the  one  given  in  Eq.  (23)  can  be 
arrived  at  for  the  coefficients  by  replacing  f  in  Eq.  (23)  by  the 
observation  vector  y  =  [t/i . . .  2/tv]T-  The  estimate  at  x  =  xo 
is  then  evaluated  to  be 

2/o  =  K°)  =  [<M°) . . .  <j>k( 0)]a 

=  [1  0  ...  0](T>TWT>)-1<f-TWy.  (25) 

Two  important  measures  of  the  quality  of  this  estimator  are 
its  bias  and  its  variance,  which  are  defined  as  follows 

•  Bias(t/0)=-E’{t/o}  -  Vo- 

•  Variance=£,{(yo  -  E{y0})2}. 

These  measures  are  controlled  by  the  two  parameters  used  in 
the  regression  procedure: 

1)  the  scale  of  the  weighting  function,  also  referred  to  as 
bandwidth  in  the  statistics  literature  (not  to  be  confused 
with  the  bandwidth  of  the  filter  G(u)  obtained  in  sub¬ 
section  II-A),  and 

2)  the  order  of  the  truncated  Taylor  series  expansion. 

Both  these  parameters  have  a  critical  effect  on  the  local 
regression  fit.  If  the  scale  is  too  small,  insufficient  data  will 
fall  within  the  chosen  support  of  the  weight  function  and  a 
large  variance,  or  a  noisy  fit  will  result.  On  the  other  hand, 
if  the  scale  is  too  large,  the  local  polynomial  may  not  fit  the 
data  well  within  the  window  of  the  weight  function,  resulting 
in  a  larger  bias.  This  is  especially  true  in  regions  where  the 


data  changes  rapidly  or  has  large  curvature.  The  order  of 
approximation  has  the  opposite  effect  on  bias  and  variance. 
A  higher  order  of  approximation  means  that  the  estimate 
fits  the  data  well,  resulting  in  a  lower  bias.  However,  since 
the  number  of  unknowns  (coefficients)  involved  in  this  is 
larger,  the  variance  of  the  estimate  is  higher.  Thus,  these  two 
parameters  must  be  chosen  to  compromise  this  bias-variance 
tradeoff. 

In  image  superresolution,  variance  is  a  measure  of  the  noise 
in  the  output  and  bias  is  a  measure  of  the  amount  of  blur 
introduced  by  the  MLS  process.  Thus,  a  technique  is  needed 
that  would  permit  selection  of  the  scale  and  order  locally  so 
that  an  effective  compromise  between  blur  introduced  by  the 
process  and  noise  in  the  output  is  reached.  This  is  really  a 
model  selection  problem  which  has  been  explored  extensively 
in  the  statistics  and  signal  processing  literature.  One  possible 
solution  is  described  in  the  next  subsection. 


B.  Selecting  Scale  And  Order.Model  Selection  Criteria 

In  order  to  evaluate  the  effect  of  scale  and  order  selection 
on  the  estimated  value,  it  is  required  to  define  a  criterion  that 
assess  the  performance  of  the  fit.  Using  mean  square  error 
(MSE)  as  a  criterion  is  not  advisable,  as  the  estimator  that 
fits  the  noisy  data  perfectly  would  then  be  the  best  estimator, 
which  is  clearly  not  desired.  The  desired  criterion  should 
indicate  poor  performance  in  case  of  both  high  bias  and  high 
variance  situations.  Two  commonly  used  criteria  are  mentioned 
below  [17]: 

1)  Prediction  mean  square  error  (PMSE)  is  defined  as 

P MSE  =  E{(ynew  -  p(xnew))2}. 

PMSE  measures  the  quality  of  prediction.  However, 
since  we  have  only  one  set  of  data,  an  estimate  of 
the  PMSE  is  implemented  by  the  cross  validation  (CV) 
estimator  [18]  which  is  defined  as 

1  N 

CV  =  —  'y  '  ( ym  —  A— m(xm))  ) 

m=  1 

where  /t_m(xm)  denotes  the  estimate  of  p(xm)  ob¬ 
tained  by  leaving  out  xm  and  computing  it  from  the 
remaining  N  —  1  points.  A  lower  value  of  CV  indicates 
a  better  fit. 

2)  The  scaled  sum  of  squared  error  (SSE)  is  defined  as 

1  N 

SSE  =  ,}  (/z(xm)  -  /)(xm))2. 

m—l 

Note  that  this  is  different  from  MSE,  which  is  ^  (t/TO  — 
/t(xm))2.  It  attempts  to  measure  how  accurately  the  esti¬ 
mate  /t(xm)  approximates  /j(xm),  instead  of  measuring 
how  well  the  estimated  points  /t(xm)  fit  the  data  ym. 
The  Cp  criterion  provides  an  unbiased  estimate  of  the 
SSE  and  is  given  by  [19] 

1  N 

Cp  =  —  ^{y-m  -  A(Xm))2  -  N  +  2p,  (26) 

m—l 
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where  p  is  the  number  of  parameters  in  the  model. 
Since  the  model  used  here  has  coefficients  ao,  ■  ■  ■ ,  a,k, 
therefore,  p  is  equal  to  k+1.  A  low  value  of  Cp  indicates 
a  better  fit.  It  should  also  be  noted  that  Cp  criterion  is 
a  modified  version  of  the  AIC  criterion  developed  by 
Akaike  several  years  back  and  used  recently  in  [20], 
For  evaluation,  the  Cv  metric  requires  much  fewer  computa¬ 
tions  than  the  CV  metric.  Thus,  in  the  simulations  and  results 
that  follow,  the  Cp  criterion  will  be  used  for  selection  of  model 
parameters.  Before  proceeding  further,  it  should  be  noted  that 
the  expression  for  Cp  given  in  Eq.  (26)  is  actually  valid  for 
“global”  regression  techniques.  For  local  regression  involving 
the  Gaussian  weight  function  used  here,  the  definition  needs 
to  be  modified  to  [21] 


Cp  = 


,  N 

1  x  -  m 

t.rCW)  a, 

v  '  m= 1 


tr(  M2) 
tr( W)  ’ 


(27) 


where, 

M2  =  ($TW$)“1$TW2$, 


and  W  is  the  diagonal  matrix  in  Eq.  (21b)  involving  the 
weights  wm  =  in(xm  —  x0).  Note  that  the  expression  for 
the  localized  Cp  requires  an  estimate  of  the  noise  variance 
cr2.  The  variance  of  the  intensity  values  in  a  region  of  the 
image  that  does  not  contain  any  sharp  edges  or  rapid  repetitive 
variations  (textures)  is  estimated.  For  regions  that  are  fairly 
constant  in  the  original  image,  this  variance  can  almost  entirely 
be  attributed  to  the  noise. 


C.  Computational  complexity 

One  of  the  advantages  of  using  local  regression  methods 
such  as  MLS  is  that  they  involve  solving  a  large  number  of 
small  system  of  equations  as  against  a  single  large  system 
usually  required  by  global  approaches.  An  analysis  of  the 
computational  complexity  of  MLS  is  provided  in  [13],  [22], 
It  is  shown  in  [22]  that  the  complexity  of  the  approximant  at 
a  single  point  x  is  bounded  by  0(Q3  +  Q'2\I(x)\  +  Q\I(x)\), 
where  Q  is  the  dimensionality  of  the  polynomial  space  and 
|/(a;)|  is  the  number  of  points  in  the  neighborhood  /( x)  of 
the  point  x  used  in  the  computation.  Clearly,  Q  and  \I(x)\ 
depend,  respectively,  on  the  order  and  scale  parameters. 

This  bound  is  valid  for  a  fixed  scale  and  order.  In  the 
adaptive  approach,  it  is  desired  precisely  to  find  a  scale  and 
order  that  would  result  in  the  optimal  tradeoff  between  blur 
and  noise.  Thus,  statistical  techniques  have  been  adopted  to 
adaptively  choose  the  scale  and  order  for  each  point.  One 
of  the  techniques  which  involves  the  Cp  metric  has  been 
described  here.  Due  to  the  complex  interdependency  between 
scale  and  order,  there  is,  so  far,  no  known  technique  that  can 
search  simultaneously  for  the  optimal  scale  and  order.  The 
usual  practice  is  to  fix  one  of  the  parameters  and  optimize 
for  the  other  and  repeat  this  process  for  different  values  of 
the  first  parameter.  This  is,  obviously,  a  more  computation 
intensive  process  than  the  simple  MLS  process  using  fixed 
scale  and  order,  since  it  involves  repeating  the  process  for 
different  values  of  scales  and  orders.  The  number  of  com¬ 
putations  performed  can  be  limited  by  restricting  the  range 


of  scale  and  order  values  over  which  the  search  is  done.  The 
justification  behind  this  is  that  real  images  can  usually  be  well 
approximated  by  piecewise  polynomials  of  orders  less  than 
five,  and  it  makes  no  sense  to  choose  the  scale  value  so  large 
that  points  very  far  away  from  the  point  being  estimated  end 
up  being  used  in  the  computation.  This  is  especially  true  in 
the  MLS  approach,  where  the  weight  function  w(x),  ideally, 
vanishes  for  arguments  x  =  X\  and  x  =  X2  with  ||a;i  —  X2W 
greater  than  a  certain  threshold  [13,  pp.35]  .There  exist  several 
other  technical  devices  to  speed  up  computations.  These  are 
based  on  the  idea  that  scale  should  be  a  smooth  function  of  the 
point  x  being  approximated  [13,  pp.43]  and  the  order  needs 
to  be  chosen  from  a  set  of  low  cardinality.  Thus,  instead  of 
evaluating  the  scale  and  order  at  all  desired  points,  these  can 
be  evaluated  for  a  smaller  set  of  points.  The  value  of  the  scale 
and  order  at  other  points  can  be  obtained  by  interpolating 
the  values  obtained  for  the  smaller  set.  Details  on  methods 
to  implement  such  techniques  efficiently  can  be  found  in  [20] 
[23], 


IV.  Simulation  Results 

Simulation  results  demonstrating  some  of  the  important 
properties  of  MLS  and  its  application  to  superresolution  are 
presented  here.  Both  synthetic  and  real  images  are  considered. 
The  PSNR  and  Msvd  metrics  are  used  as  indicators  of 
performance.  The  Msvd  metric,  defined  in  [5],  is  based  on 
the  singular  values  of  the  image  and  gives  a  much  better 
representation  of  the  quality  of  an  image  with  regard  to 
the  human  visual  system  for  a  large  variety  of  distortions. 
The  Msvd  metric  has  also  been  used  in  the  context  of 
second-generation  wavelet  superresolution  [4]  [6]  to  choose 
a  threshold  value  for  the  wavelet  coefficients  that  strikes  an 
optimal  balance  between  noise  filtering  and  blur  introduced 
by  the  thresholding.  A  lower  value  of  Msvd  indicates  better 
image  visual  quality. 

Figure  3  shows  four  LR  frames  of  size  64  x  64  of  the 
synthetic  “cameraman”  image  corrupted  by  noise  of  variance 
0.01.  The  desired  resolution  enhancement  factor  is  2;  hence  the 
resulting  HR  images  are  128  x  128.  The  HR  images  obtained 
by  choosing  a  fixed  scale  and  order  for  the  entire  image, 
along  with  their  corresponding  PSNR  and  Msvd  metrics  are 
shown  for  different  values  of  scales  and  orders  in  Figures  4 
and  5.  The  tradeoff  between  blur  (bias)  and  noise  (variance) 
described  in  subsection  II-C  can  be  clearly  seen  from  these 
images.  The  HR  image  obtained  by  adaptively  selecting  the 
scale  and  order  is  shown  in  Figure  6.  Even  though  the  image 
obtained  by  the  adaptive  procedure  does  not  have  the  highest 
PSNR,  or  the  lowest  Msvd ,  certain  distinctive  features  in  the 
image  such  as  the  leg  of  the  tripod  and  the  building  in  the 
background  are  captured  much  better  in  the  adaptive  mode  of 
image  reconstruction  than  in  any  of  the  fixed  scale  and  order 
modes. 

Next,  the  results  obtained  by  using  MLS  in  superresolution 
are  compared  to  those  obtained  by  the  surface  approximation 
method  based  on  Delaunay  triangulation  [7].  Figure  7(a)  shows 
one  sample  from  the  six  64  x  64  size  LR  frames  of  the  tank 
image  sequence,  corrupted  by  noise  of  variance  0.005.  Here 
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(a) 


(b) 


(0 


Fig.  4.  Reconstructed  images  with  approximation  orders  0,2,4  for  a  fixed  value  of  scale;  (a)  Scale=0.4,  Order=0,  PSNR=22.78,  MgvD=^- 08;  (b)  Scale=0.4, 
Order=2,  PSNR=23.07,  MSVd=  0.05;  (c)  Scale=0.4,  Order=4,  PSNR=20.25,  MSVd=  0.16 


(a) 


(b) 


(c) 


Fig.  5.  Reconstructed  images  with  different  chosen  scales  for  a  fixed  order  of  approximation;  (a)  Scale=0.4,  Order=4,  PSNR=20.25,  M$vd=  0.16;  (b) 
Scale=0.6,  Order=4,  PSNR=22.93,  MSvd=  0.06;  (c)  Scale=0.8,  Order=4,  PSNR=23.53,Msyr>=0.05 


again,  the  resolution  enhancement  factor  is  2.  The  resulting 
HR  image  obtained  from  the  Delaunay  triangulation  method 
is  shown  in  Figure  7(b)  and  that  obtained  from  the  MLS 
technique  is  shown  in  7(c).  It  is  seen  that  owing  to  the  adaptive 
nature  of  MLS,  it  performs  better  in  filtering  noise  while  at 
the  same  time  preserving  edges. 

Finally,  the  result  of  applying  the  Moving  Least  Squares 
method  to  a  real  data  sequence  (supplied  by  Air  Force  Re¬ 
search  Laboratory  at  Rome,  NY)  is  shown  next.  Figure  8(a) 
shows  one  of  the  six  90  x  90  size  LR  frames  obtained  from 
a  video  sequence.  Note  that  since  this  is  a  real  data  sequence 
captured  by  a  video  camera,  it  has  some  amount  of  input  blur. 
To  perform  registration,  it  is  necessary  to  know  the  motion 
parameters  of  the  camera.  Since  these  are  not  known  a  priori, 
these  are  estimated  by  the  projective  model  in  [14].  The  results 
of  applying  the  Delaunay  triangulation  method  and  the  MLS 
technique  are  shown,  respectively,  in  Figures  8(b)  and  8(c). 
The  resulting  HR  images  obtained  by  both  methods  are  seen 
to  be  comparable  in  visual  quality.  Also,  the  HR  images  from 
both  methods  are  seen  to  contain  a  small  amount  of  blur.  This 
is  to  be  expected  since  no  post  deblurring  has  been  performed 


Fig.  6.  Reconstructed  image  by  adaptive  selection  of  scale  and  order. 
PSNR=22.50,Msv^=0.06 

on  the  HR  image. 

V.  Conclusions 

In  the  continuous  formulation  of  MLS  applied  to  noise 
filtering  and  approximation  of  irregularly  spaced  data  in  the 
superresolution  problem,  an  expression  for  filter  bandwidth  in 
terms  of  the  order  of  polynomial  approximation  and  the  scale 
of  the  Gaussian  weighting  function  is  obtained.  The  discrete 
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Fig.  7.  Simulation  results  on  tank  sequence  showing  comparison  between  Delaunay  Triangulation  and  MLS  methods,  (a)  Sample  LR  frame  of  the  tank 
sequence;  (b)  HR  image  using  Delaunay  Triangulation.  PSNR=25.02,  Msvd- 0.04;  (c)  HR  image  using  MLS.  PSNR=31.22,  Msvd=Q-QQ7  ■ 


(a)  (b)  (c) 

Fig.  8.  Simulation  results  showing  comparison  between  Delaunay  Triangulation  and  MLS  methods  on  a  real-world  data  sequence,  (a)  Sample  LR  frame  of 
the  building  sequence;  (b)  HR  image  using  Delaunay  Triangulation;  (c)  HR  image  using  MLS. 


implementation  of  MLS  is  performed  on  a  sequence  of  LR 
images  registered  to  form  the  grid  of  non-uniformly  spaced 
data  points.  The  MLS  technique  is  then  used  adaptively  to 
evaluate  the  values  of  the  high-resolution  image  on  a  regularly 
spaced  high-resolution  grid.  The  parameters  for  adaptation  in 
the  local  approximation  are  the  scale  and  order.  Noise  in  the 
output  image  increases  with  increasing  order  of  approximation 
but  reduces  with  increasing  scale.  One  the  other  hand,  the  blur 
introduced  reduces  with  increasing  order  but  increases  with 
increasing  scale.  For  a  given  scale  of  the  Gaussian  weight 
function,  the  PSNR  first  increases  with  order  of  the  image 
and  then  decreases  again.  This  is  because  at  low  orders, 
degradation  due  to  blur  is  higher  than  that  due  to  noise, 
while  at  higher  orders, the  reverse  happens.  Between  these 
two  extremes,  a  tradeoff  between  the  degrading  effects  of 
blur  and  noise  produces  a  better  image.  This  suggests  that 
there  is  an  optimal  value  of  approximation  order  for  a  given 
scale.  Similarly,  there  is  an  optimal  value  of  scale  for  a  given 
approximation  order.  Statistical  and  signal  processing  criteria 
for  joint  optimization  of  the  two  dependent  parameters  are 
applied  to  get  optimal  local  tradeoff  between  noise  filtration 
and  reduction  of  blur  introduced  by  the  MLS  process. 

The  availability  of  two  parameters  has  the  capability  of 
producing  greater  robustness  in  numerical  implementation  than 
when  only  one  parameter  (threshold)  is  used  to  achieve  optimal 


tradeoff  between  noise  filtering  and  blur  removal  [4],  Ideally, 
it  would  be  desirable  to  introduce  estimate  of  blur  type  and 
blur  support  parameters  [7]  [8,  and  references  therein]  into  the 
superresolution  process  with  not  only  noise  filtering  but  also 
simultaneous  deblurring. 


One  goal  of  future  research  is  to  study  the  approach 
advanced  when  input  blur  is  present.  Special  attention  may  be 
required  to  obtain  the  irregularly  sampled  raster  in  Figure  1  in 
the  presence  of  near  field  and  medium  field  blurs  for  general 
camera  motion.  The  second  goal  is  to  study  the  problem  in 
the  setting  of  multivariate  orthogonal  polynomials  that  are 
not  product  separable  and  by  using  weight  functions  that  do 
not  necessarily  have  circular  support.  For  example,  if  the  n- 
dimensional  weight  function  is  chosen  to  be  the  multivariate 
Gaussian  given  by 


w(x)  = 


1 

(27r)n/2<Ti  .  .  .  On  6  P 


1  (X1 

2  W 


then  a  set  of  (/-dimensional  non-separable  orthogonal  Hermite 
polynomials  can  be  found,  as  was  proved  by  Grad  in  [24] 
for  the  case  when  ui  =  02  =  •  •  •  =  crn  =  1.  The  filter 
impulse  response  in  the  n-dimensional  case  can  be  obtained 
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by  generalizing  the  expression  in  Eq.  (30)  as, 

fif(x)  =  w(x)Y^  (-j)  \\H^  (x;  <T) ' 

i=0  '  ' 

where, 

X  =  [ari  x2  ■  ■  ■  xn]T,  i  =  [*1  *2  . . .  in]T, 

|i|  =  *1  +  *2  +  •  •  •  +  i-  = 

1  l\  I2  In 

E=EE-E 

i— O  2i=0  42—0  iri— 0 

and  -ffJ,;1'  (x;  <x)  is  the  multivariate  non-separable  Hermite 
polynomial  of  total  degree  |2i|  with  the  degree  of  x:l  given 
by  ij  for  j  =  1,2, ...  ,  n.  Adapting  the  shape  of  the  weight 
function,  and  hence  the  neighborhood  selected  for  each  point, 
is  shown  to  be  a  very  effective  technique  in  reducing  blur 
[25],  However,  the  technique  described  in  [25]  does  not  adapt 
the  order  of  approximation  at  each  point.  It  is  expected  that 
by  including  such  shape  adaptivity  with  the  scale  and  order 
adaptivity  presented  here,  much  better  results  can  be  obtained. 
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Appendix  I 

Derivation  of  the  filter  bandwidth 


A  derivation  of  the  bandwidth  of  the  filter  whose  impulse 
response  is  given  by  Eq.  (7)  is  presented  here.  For  odd  i, 
Hi(x)  is  an  odd  function  of  x  and  hence  </>i(0)  =  0  .  As  a 
result,  only  those  terms  in  the  summation  in  Eq.  (7)  need  to 
be  retained  for  which  i  is  even.  The  summation  in  (7)  can, 
therefore,  be  rewritten  as 


sO)  = 

2= 0 


(t>2i{$)(t>2i{x) 

22i(2  z)! 


w{x), 


(28) 


where  l  =  k/ 2,  if  k  is  even  and  (  =  (k  —  l)/2  if  k  is  odd.  It 
can  be  shown  that, 

J?2n(0)  =  (-1)"^. 

Using  this,  and  the  expressions  for  w(x)  and  4>i(x),  Eq.  (28) 
simplifies  to 


a(x)  =  2 L  — 4<ii — «>(*) 


i= 0 


e-2/2-2  ' 


(29) 

(30) 


Interestingly,  the  preceding  form  for  g(x)  can  be  extracted 
from  a  result  in  [26],  which  is  derived  by  considering  the  effect 
of  the  quantum  free  propagator  on  wave  packets  approximated 
by  a  linear  combination  of  Hermite  polynomials. 


Replacing  x  by  x/cr\/2  in  Eq.  (4)  and  using  the  result  in 
(29)  yields  after  simplification. 


9(x)  = 


2v(~l)i  d2i 


i= 0 


2 li\  dx 2l 


OO)) . 


The  FT  G(u)  of  g(x)  is 


2—0 


2  —  0 


Therefore, 


1  1 


J-00  i=0  j=o  J-°°  J ' 

(31) 

The  FT  W(oj)  of  tv(x)  is  easily  shown  to  be  W( ui)  = 
e-w2CT-’/2  Substituting  this  in  Eq.  (31), 


\  1 2 j  4-^  r  (*2“2)i+j 

|G(a;)|du;  =  ^W  e 

i=0j=oJ-°°  J ' 

'  '  rrW+ni 

2—0  j  —  0 

EEj::.;;:1',,.  (32) 


du 


>[2  (i  +  j)]'. 


V7T 

a 


2  —  0  j—0 


2*+ii!j!4*+j0-2(i+i)  (j  _|_  j)i 

[2  (i  +  M 
8  i+H\j\(i  +  j)V 


since 


=  1.3  .  .  .  (2n  -  1)  /  7T 


2n<72 


(2  n)\ 


a 2  er2ra4”n!  V  cr2 
(33) 


Similarly, 


w 


■\G(u)\2duj  - 


(2i  +  2j  +  l)[2(i  +  j)]\ 


2 a3  8 l+H\j\(i  +  j)! 


(34) 


Substituting  Eq.  (32)  and  Eq.  (34)  in  the  expression  for 
bandwidth. 


A  ,  = 


njGMi2^ 


one  gets. 


A2  ,= 


\fn  (2i+2j  +  l)[2(i+j)]! 

2cr3  Z-/i=0  Z^j=0  8t+J’i!j!(i+j)! 


1 

2a2 


■y/TT  W 

a  Z-/2=0  2-/  j=i 

1  + 


Mi+j)}'. 
i= 0  Z—*j=0  i!j!(i+j)! 

y*  42(»+l)]! 

2-^i=0  0  8i+i i\j\(i+j)\ 

[2(i+j)V- 

0  0  8i+H\j\(i+j)\ 


El  \r~~\  l 
2—0  2-/ 7=1 


9  j[2(2+j)]! 

z  Z^2=0  7=0 


=0  Z-/j= 0  8*+^2!j!(2+j)! 


El 

2  —  0  2^7  = 


[2(»+j)]! 


1 

2^2 


0  Z-^j=0  8i+-?2!j!(2+j)! 
I 


i  + 


iV  y;  IMi+Hll 

Z-/2=0  Z-/j= 0  8i+Ji!j!(2+j)! 


El  sr~~^  l 

2—0  2-/  7=i 


[2(»+j)]t 


2=0  Z-/j=0  8*+^2!j!(2+j)! 


(35) 
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We  believe  that  the  formula  for  Af,  derived  above  is  new  and 

provides  an  analytical  basis  for  guidelines  in  the  parameter 

selection  process. 
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